J I Cardesa
email: jcardesa@onera.fr

R Fiévet

E Piot

H Deniau

C Airiau

Optimizing an acoustic liner by automatic differentiation of a compressible flow solver

Keywords: gradient-based multi-parameter optimisation, automatic differentiation, time-domain impedance boundary condition, acoustic liner

The attenuation power of an acoustic liner is optimized by unsteady compressible flow simulations that model the presence of the liner with a time-domain impedance boundary condition (TDIBC). Two test cases are considered: a one-dimensional impedance tube and a two-dimensional grazing incidence tube. The impact of the liner on the pressure field is tuned through various TDIBC parameters imposed at the boundary where the liner surface is located. Automatic differentiation of the direct numerical simulation (DNS) code yields the gradients of the time-integrated pressure perturbations with respect to the TDIBC parameters. The gradients are subsequently used in an iterative process that minimizes the perturbations at a given location. In both test cases considered, TDIBC parameters are found which are more effective at attenuating perturbation frequencies other than those the initial TDIBC is designed to cancel. Finally, a sensitivity analysis reveals the importance of the choice of the initial TDIBC parameters for finding optimized settings. Our results pave the way to a wider application of the presented procedure in more complex problems with bulk flow motion and larger perturbation amplitudes.

Introduction

Optimizing liners in order to improve their attenuation capacity is the subject of much experimental and theoretical work. One strategy that receives significant attention is that of segmented liners. In early numerical studies, Baumeister [START_REF] Baumeister | Evaluation of optimized multisectioned acoustic liners[END_REF] analyzed a number of axially segmented liner configurations for liners in a rectangular duct without any mean flow, and concluded that their use fails to offer sufficient advantage over a uniform liner except in low-frequency, single-mode applications. More recently, Law et al. [START_REF] Law | Optimisation of axially segmented liners for aeroengine broadband noise[END_REF] reach similar conclusions. By choosing more elaborate liner arrangement combinations which scatter sound towards different modes that are then attenuated by a subsequent stage of liners, Gerhold et al. [START_REF] Gerhold | Segmented liner to control mode scattering[END_REF] report acoustic performance improvements. Following a rather theoretical approach, Campos et al. [START_REF] Campos | On the optimization of non-uniform acoustic liners on annular nozzles[END_REF] consider wave convection within a uniform axial flow, where a circumferentially non-uniform but axially uniform wall impedance distribution allows them to speed up axial decay of the slowest decaying wave modes by tuning the circumferential impedance distribution. A different approach found in studies by Nark et al. [START_REF] Nark | Broadband liner optimization for the source diagnostic test fan[END_REF][START_REF] Nark | Modeling of broadband liners applied to the advanced noise control fan[END_REF] consists of setting a cost function based on in-duct attenuation, and finding its maximum value on a tiled reactance-resistance grid. The process is repeated at several frequencies, and a close-in step on the attenuation maximum is carried out using a finer reactance-resistance grid. The numerous noise propagation calculations induced by this space-sampling approach can be alleviated when using an optimization package external to the flow solver, as done in the more recent studies [START_REF] Nark | Further development and assessment of a broadband liner optimization process[END_REF][START_REF] Sutliff | Design and acoustic efficacy of a broadband liner for the inlet of the DGEN aero-propulsion research turbofan[END_REF]. It is worth noting that these studies rely on the Ffowcs-Williams and Hawkings analogy for the sound propagation computations. Özkaya et al. [START_REF] Özkaya | Development of an adjoint caa solver for design optimization of acoustic liners[END_REF][START_REF] Özkaya | A two-level approach for design optimization of acoustic liners[END_REF][START_REF] Özkaya | Efficient design optimization of acoustic liners for engine noise reduction[END_REF] opt for a rather different computational approach, which deserves mention since it is a direct precursor to our study. They solve the steady flow equations using a commercial computational fluid dynamics (CFD) software, and the obtained mean velocity and pressure fields are given as external input to the computational aeroacoustics (CAA) solver based on the time-dependent linear Euler equations (LEE). The acoustic liner within the CAA framework is modelled according to a 5-parameter extended Helmholtz resonator (EHR), and derivatives of a cost function based on acoustic wave propagation over a time interval are obtained by automatic differentiation (AD) of the CAA solver. Finally, following a gradient descent algorithm, an optimal set of liner parameters can be found. Such an approach has the distinct advantage of allowing for an adjoint-based AD optimization in space that yieds the high-dimensional design gradient vector in one stroke.

The study presented here contrasts with the work on liner optimization previously mentioned in that a fully non-linear and unsteady Euler/Navier-Stokes solver is used without decoupling the flow from the noise computations. Our approach is based on the solution of the governing equations without resorting to acoustic analogies or the linearized Euler equations. Hence, AD is used directly on the DNS flow solver. What is similar to the work of Özkaya et al. [START_REF] Özkaya | Development of an adjoint caa solver for design optimization of acoustic liners[END_REF][START_REF] Özkaya | A two-level approach for design optimization of acoustic liners[END_REF][START_REF] Özkaya | Efficient design optimization of acoustic liners for engine noise reduction[END_REF] is that derivatives of a cost function with respect to the acoustic liner parameters are computed, with the cost function to be minimized being defined over a time interval. Unlike Özkaya et al. [START_REF] Özkaya | Development of an adjoint caa solver for design optimization of acoustic liners[END_REF][START_REF] Özkaya | A two-level approach for design optimization of acoustic liners[END_REF][START_REF] Özkaya | Efficient design optimization of acoustic liners for engine noise reduction[END_REF], the flow solver is unsteady, relies on a high-order spatial discretization, and the liner is not modelled on the basis of the EHR model. Instead, the oscillo-diffusive representation [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] of the reflection operator at the liner surface is used to model the liner through a time-domain impedance boundary condition [START_REF] Fiévet | Strong compact formalism for characteristic boundary conditions with discontinuous spectral methods[END_REF] (TDIBC). As a counterpoint, it must be stated that -at this stage -the problems considered are simpler in terms of geometric complexity than those treated by Özkaya et al. [START_REF] Özkaya | Development of an adjoint caa solver for design optimization of acoustic liners[END_REF][START_REF] Özkaya | A two-level approach for design optimization of acoustic liners[END_REF][START_REF] Özkaya | Efficient design optimization of acoustic liners for engine noise reduction[END_REF], and limited to a design gradient vector of much lower dimensionality -no spatial variability of the liner properties is considered here.

It is important to note that the cases analyzed in this study are chosen so as to place the tools that have been developed in a validation and proof-of-concept environment. The problems that are solved involve small pressure perturbations in the absence of bulk flow, and hence LEE should suffice to obtain the acoustic field. But the methods employedgradient-descent optimization through AD of a non-linear viscous solver -are indifferent to the physical problem considered. It follows that acoustic propagation cases where the LEE assumptions break down are still within reach of our approach. Such flow problems 2 are left for future work, and the present study intends to describe the computational procedure and assess its performance in controlled scenarios where it is expected that improved liner properties can be obtained. The computational tools will be described first in Section 2, followed by the test cases in Section 3. Results will be outlined in Section 4 and discussed in Section 5, before conclusions are brought forward in Section 6.

Computational Tools

Numerical flow solver

The flow solver used in the present study is called JAGUAR [START_REF] Cassagne | Jaguar: a new CFD code dedicated to massively parallel high-order LES computations on complex geometry[END_REF][START_REF]Comparison of various CFD codes for LES simulations of turbomachinery: from inviscid vortex convection to multi-stage compressor[END_REF], which is jointly owned and developed by CERFACS1 and ONERA. This CFD code uses a spectral difference scheme [START_REF] Liu | Spectral difference method for unstructured grids I: basic formulation[END_REF][START_REF] Kopriva | A staggered-grid multidomain spectral method for the compressible navier-stokes equations[END_REF][START_REF] Jameson | A proof of the stability of the spectral difference method for all orders of accuracy[END_REF][START_REF] Van Den Abeele | On the stability and accuracy of the spectral difference method[END_REF][START_REF] Vanharen | Revisiting the spectral analysis for high-order spectral discontinuous methods[END_REF] for the spatial discretization, which belongs to the family of finite element methods. It solves the strong form of the governing equations for a compressible flow. The method allows for an arbitrarily high order of accuracy even on unstructured meshes, and is suited to both compressible or incompressible flow computations. In practice, the scheme reconstructs a polynomial solution of the conservative variables U = (ρ, ρu, ρv, ρe) ⊺ inside each element with density ρ, streamwise (herein x-coordinate) velocity u, wall-normal (herein y-coordinate) velocity v and specific internal energy e. The polynomial order p determines the overall scheme's accuracy and computational effectiveness. Throughout this study, p = 4, which results in a fifth-order accurate spectral difference scheme. Furthermore, the time integration follows an optimized low-dissipation low-dispersion fourth order 6-step Runge-Kutta scheme inspired from Berland et al. [START_REF] Berland | Low-dissipation and low-dispersion fourth-order Runge-Kutta algorithm[END_REF]. Finally, all results presented in the following sections are carried out in inviscid (Euler flow) mode in one or two dimensions, so that the governing equations in generalized coordinates (ξ,η) are given by:

1 J ∂U ∂t + ∂ ∂ξ E c ξ x J + F c ξ y J + ∂ ∂η E c η x J + F c η y J = 0, ( 1
)
where J is the mesh local Jacobian and the convective fluxes of U in physical space are called (E c , F c ). We call p the static pressure such that the fluxes read as:

E c = (ρu, ρu 2 + p, ρvu, (ρe + p)u) ⊺ (2a)
F c = (ρv, ρuv, ρv 2 + p, (ρe + p)v) ⊺ (2b)
Choosing to set JAGUAR to its inviscid rather than its viscous mode is justified on the basis of the test cases considered which lack bulk flow motion. However, the procedure of automatically differentiating JAGUAR is indifferent to the code being set to viscous or inviscid mode, so that a broader set of flow problems can be considered in the future.

Time-domain impedance boundary conditions

In order to capture the effect of the acoustic liner on the flow field, a TDIBC is used as developed in Fiévet et al. [START_REF] Fiévet | Strong compact formalism for characteristic boundary conditions with discontinuous spectral methods[END_REF]. It consists of three separate components detailed below as a bulleted list2 :

• An impedance law z(s), which characterizes the acoustic properties of the coating in Laplace-space (s = jω, where j is the complex unit number and ω the angular frequency). For practical reasons, we work with the reflection coefficient B defined in Laplace space as:

B(s) = z -ρc z + ρc ( 3
)
where c is the speed of sound. Importantly, the modulus of B is always inferior or equal to unity, which offers interesting stability properties to a B-based TDIBC as opposed to a z-based one, as shown by Monteghetti et al. [START_REF] Monteghetti | Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized euler equations[END_REF].

• A discretized representation of B, referred to as β. This representation consists in a delayed multi-pole approximation [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] which models a single-degree of freedom acoustic liner -the kind used throughout the present study. It relies on N s oscillatory and N ν diffusive poles, and is defined as follows:

β(s) = β ∞ + Ns n=1 r 1,n s -s n + Nν k=1 µ 1,k s + ν k + e -sτ Ns n=1 r 2,n s -s n + Nν k=1 µ 2,k s + ν k . (4)
β depends on a parameter vector q that includes the bulk reflectivity β ∞ , a time delay τ related to the depth of the liner cavities, the oscillatory weights r 1,n and r 2,n , the oscillatory poles s n , the diffusive weights µ 1,k and µ 2,k , and the diffusive poles ν k . The more poles are used, the more broadband-accurate the representation can become.

• A semi-discrete formulation of this impedance condition in the time-domain to couple itself with the internal flow solver. In JAGUAR, the TDIBC is implemented in characteristic form using a time-domain scattering operator B. This amounts to evaluating the amplitude of the incoming characteristic at the boundary as a function of the outgoing characteristic according to:

∂ t p + (ρc)∂ t u.n = B[∂ t p -(ρc)∂ t u.n] = β ⋆ [∂ t p -(ρc)∂ t u.n], (5)
where u.n is the velocity normal to the boundary and β is the inverse Laplace transform of β. The operator B is evaluated by transforming the convolution product into a sum of (N s +N ν ) auxiliary variables, i.e. one per acoustic pole, each solution of an ordinary differential equation [START_REF] Monteghetti | Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized euler equations[END_REF] defined by the parameters inside vector q.

A thorough presentation of the TDIBC development is given in Monteghetti et al. [START_REF] Monteghetti | Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized euler equations[END_REF]. For a validation and verification of the numerical implementation of the TDIBC, the reader is referred to Fiévet et al. [START_REF] Fiévet | Strong compact formalism for characteristic boundary conditions with discontinuous spectral methods[END_REF], which includes a detailed explanation of its deployment within a spectral difference flow solver such as JAGUAR. With regards to the present study, the important point is that the TDIBC is an explicit function of a finite set of both real and complex parameters contained inside q.

Automatic differentiation

Gradient-based methods play an important role in the process of seeking those conditions that modify the solution of a physical problem towards a specific goal, either for optimization or control purposes. Computing the necessary gradients (also referred to as sensitivities) is a task that can be automated by automatic differentiation (AD) tools. In short, AD augments a given "primal" code that initially computes outputs Y i from inputs X j into a "differentiated" code that additionally computes some derivatives dY i /dX j requested by the user. AD provides two main modes, the tangent/direct/forward mode and the adjoint/reverse/backward mode. If i and j are the indices along the output and input spaces, respectively, bound so that 1 ≤ i ≤ m and 1 ≤ j ≤ n, then the tangent mode is most efficient when n ≪ m while the adjoint mode is the only realistic option for m ≪ n.

Over the years, the AD community has produced a number of tools to generate differentiated code based on different strategies, chiefly operator overloading (OO) and source transformation (ST). Perhaps the most important point to bear in mind when choosing an AD tool is the programming language in which the primal code is written. JAGUAR is written following features from the Fortran 90 standard onward, which narrows down our choice of an AD tool towards an effective and well established ST-based open-source software: TAPENADE [START_REF] Hascoët | The Tapenade automatic differentiation tool: principles, model, and specification[END_REF]. In a previous study [START_REF] Cardesa | Adjoint computations by algorithmic differentiation of a parallel solver for time-dependent pdes[END_REF], the compatibility between TAPE-NADE and JAGUAR was verified and the performance of the differentiated code in both tangent and adjoint modes was assessed for gradient computations. In particular, it was shown that the computational cost of the adjoint-differentiated code was above 15 times that of the primal code, while for the tangent-differentiated code the same factor was reduced to 1.7. Both numbers refer to the cost of running JAGUAR in explicit time-integration mode to compute one sensitivity. The conclusion from those numbers is that the adjoint-mode differentiation begins to pay off when more than 15 sensitivities are required in a single code execution. Since the flow problems that are solved further down involve computing gradients in 14-dimensional space at most, the (vector) tangent mode is chosen. For ease of implementation, the serial version of JAGUAR is chosen for differentiation over the parallel version.

Problem Setup

In the present section, the two considered test cases are described and the common optimization procedure is explained. Both test cases are computed assuming inviscid flow, time-integrated using a constant time step size, and spatially-discretized with a fifth-order approximation (p = 4).

Wave packet in a 1D impedance tube

A one-dimensional problem representing a tube of length L = 1.8 m is considered, open at x = 0 and closed with a TDIBC at x = L. A sketch is provided in Fig. 1. The pressure in the tube is at its ambient value p a = 101325 Pa. At t = 0, a pressure wave packet enters the tube at x = 0, travelling at the speed of sound c = 344.37 m/s towards x = L. The properties of the wave packet are illustrated in Fig. 2 and detailed in its caption. It is emphasized that the input pressure perturbations are broadband throughout the present study. The spatial discretization relies on N x = 60 grid points, 2 . For all cases considered in the present study, A = 1 Pa, σ 2 = 400. The two figures above were based on f = 1000 Hz, but input wave packets with f = 1500 Hz and f = 2000 Hz are also employed in this study.

′ = A sin (2πf (t -3/f )) exp -σ 2 f (t -3/f )
resulting in a constant grid spacing of ∆x = 3.051 × 10 -2 m and a Nyquist frequency of c/(2∆x) = 5643 Hz. As can be seen in Fig. 2(b), the amplitude spectrum of the input wave packet shows no energy well before reaching the Nyquist frequency.

A probe at x = 0.3 m records the pressure fluctuations p ′ with respect to p a observed in Fig. 3. The pressure fluctuations induced by the input wave packet are those recorded during the interval between t = 0 and t 1 = 8 × 10 -3 s. The second wave packet recorded by the probe during t 1 ≤ t ≤ t 2 corresponds to the pressure waves reflected by the TDIBC while they travel back towards x = 0. The simulation ends at t 2 , reached by setting a constant time step ∆t = 7 × 10 -7 s and running for 25000 time steps. It is worth emphasizing that the ability to resolve accurately pressure perturbations of the order of 1 Pa with respect to ambient pressure, as done in the present study, illustrates the suitability of the high-order code JAGUAR at tackling acoustic propagation problems.

Wave packet grazing over an acoustic liner

The second test case is inspired on the grazing incidence tube (GIT) experiment of Jones et al. [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF]. In that experiment, a three-dimensional tube with square cross-section is used to collect reference data for validation of acoustic propagation codes.

In the present study, the domain is discretized with a two-dimensional mesh. A sketch of the setup for this test case is provided in Fig. 4. As can be seen in that figure, the presence of an acoustic liner at the bottom wall is modelled with a TDIBC that interacts with the input right-running pressure waves. A regular grid with constant ∆y = 4.5 mm and ∆x = 4.2 mm was chosen, leading to a Nyquist frequency of 37875 Hz. The same grid is used in a previous study [START_REF] Fiévet | Strong compact formalism for characteristic boundary conditions with discontinuous spectral methods[END_REF], where a harmonic plane pressure wave at 2000 Hz is input at the inlet of the domain. In the present study, a Gaussian plane wave packet is used with properties as those displayed in Fig. 2, but with a peak amplitude spectrum centered at 1500 Hz instead of 1000 Hz as that shown in Fig. 2 of the pressure fluctuations in the domain can be seen on the sequence in Fig. 12. All computations for this test case are carried out with a constant time step of 2 × 10 -6 s for 5005 time steps.

(b). The transient evolution

Optimization process

In both test cases previsouly described, JAGUAR imposes a TDIBC at a boundary by relying on a vector q of values for the coefficients of the ODR of the reflection operator found in Eq. ( 4). The initial inputs of q, which are displayed in the left-hand column of Table 1 as q ref , have been chosen to match the acoustic properties of the ceramic tubular liner (CT57) found in Jones et al. [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF] by methods outlined in the work of Monteghetti [START_REF] Monteghetti | Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized euler equations[END_REF]. When comparing Eq. ( 4) and Table 1, it can be observed that the sums in Eq. ( 4) are limited to N ν = 0 and N s = 4, reducing the number of input parameters controlling the TDIBC to 14 -the dimension of q.

For better compatibility with the automatic differentiation step, JAGUAR has to be mildly modified so that a top-level subroutine explicitly takes the 14 control parameters of vector q as input, computes the flow solution by integrating the Euler equations in time, and then returns a scalar cost function as output. The cost function is defined as follows:

J(q) = t2 t1 [p ′ (q, t) -p ′ tar (t)] 2 dt, (6)
where p ′ tar (t) is the target perturbation pressure signal, while p ′ (q, t) is the perturbation pressure signal obtained with input q. Both pressure signals are measured at a fixed location. For the 1D impedance tube, x = 0.3 m is chosen, while for the 2D GIT case the probe location is on the wall opposite to the TDIBC and between the streamwise position of the TDIBC end and the outlet: x = 70 cm and y = 5 cm -see Fig. 4. In general, optimizing the pressure fluctuations at a single point may lead to generating locally a pressure node. However, that possibility is ruled out in the flow-less test cases considered due to the fact that the outlet is perfectly non-reflecting, so that no standing waves are created after the plane waves employed as input have passed by.

The limits of the time integral in Eq. ( 6) are chosen so as to contain the temporal window during which the wave packet goes through the probe location. For the 1D impedance tube, t 1 and t 2 are as displayed in Fig. 3 4). The missing poles and their corresponding weights are the complex conjugate of already-given poles and weights: s3 = s1 , s4 = s2 , r1,3 = r1,1 , r1,4 = r1,2 , r2,3 = r2,1 , r2,4 = r2,2 . The left-hand column displays the values found in Monteghetti [START_REF] Monteghetti | Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized euler equations[END_REF] to match the experimentally-measured reactance properties in Jones [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF]. The middle column contains the values after optimization in the 1D impedance tube when all elements of q can be changed except for τ , while the right-hand column is obtained by keeping both τ and β∞ constant.

the pressure perturbations damped by the TDIBC. In the optimization procedure that follows, the goal is to find the set of input parameters q opt that brings p ′ (q, t) closest to p ′ tar , in other words q opt should minimize J. Since we aim to minimize the pressure perturbations, the chosen p ′ tar appears to be an ideal target. Of course one has no guarantee that it is realizable, yet exploring how close one may come to it is of great interest.

Once JAGUAR contains the top-level subroutine where the dependence between the input vector q and the scalar output J(q) has become explicit, the subroutine can be automatically differentiated with TAPENADE using vector tangent mode. The differentiated top-level subroutine not only outputs the scalar J(q), but also a vector dJ which contains the sensitivities of J with respect to each element in q. A validation of the sensitivities obtained with AD is performed for the 1D impedance tube by comparing them against finite-difference (FD) approximations. Both estimates of the same sensitivities are displayed in Table 2 and are found to agree reasonably well. The FD estimates are obtained by successive additional runs of JAGUAR (in its non-differentiated form) while changing one element of q at a time. Each element of q is increased by multiplying it times (1 + 10 -5 ). The agreement between FD and AD is sufficient to validate the differentiation procedure, yet perceptible differences are to be expected since FDs suffer from truncation error while AD estimates do not. In terms of CPU time, FD requires 14 additional code executions to compute 14 derivatives. The differentiated version of Jaguar was found to require the equivalent of 11.6 additional code executions to provide all 14 sensitivities -1 code execution takes approximately 56 s on a desktop computer. So not only is the AD more accurate at computing sensitivities, it is also faster -two distinct advantages when it comes to converging towards a minimum through gradient descent algorithms.

The AD-obtained sensitivities are embedded into the following steepest descent loop:

q n+1 (m) = q n (m) -λ [dJ(m)] n , ( 7
)
where m is the index of the element within vectors q and dJ, while n is the iteration number. The parameter λ is a fixed "small" number, chosen based on the rate of convergence for each test case.

The values obtained during updates of q in Eq. ( 7) need to be checked at every iteration so that they comply with the following conditions:

• the real part of the oscillatory poles s1→4 has to remain negative

• the modulus | β(s)| has to remain below or equal to unity for all explored frequencies.

In practice, we checked all frequencies between zero and 5650 Hz.

The upper bound of the frequencies checked is chosen based on the Nyquist frequency being 5643 Hz for the 1D impedance tube, and the fact that the amplitude spectrum of the input pressure perturbations is as shown in Fig. 2(b) -negligible energy is found above 3000 Hz. The same frequency check is used in the 2D GIT case since similar input wave packets were used -with f = 1500 Hz instead of 1000 Hz as in Fig. 2(b).

Results

The results obtained during the optimization under different conditions are presented next, listing simply the findings. A discussion of the results including their interpretation is delayed until Section 5.

1D impedance tube

The first problem to be tested by means of the optimization procedure outlined in subection 3.3 is the 1D impedance tube. The descent algorithm in 14-dimensional space following Eq. ( 7) converges very quickly towards a trivial solution that relies on increasing τ . Such a tendency towards attenuating the energy of the reflected pressure perturbations at x = L by lengthening the delay in the TDIBC response τ is clear from Table 2, where dJ/dτ is found to be orders of magnitude larger than the other sensitivities. However, it corresponds to an unphysical solution consisting of an infinite-depth liner. Indeed, the delay τ describes the back-and-forth time of an acoustic wave travelling at the speed of sound inside the porous material. Hence τ must remain small enough to ensure that the coating remains thin enough for practical reasons, and so it is decided to exclude τ from the tunable parameters thereby removing it from the optimization procedure. That is the reason why q ref (2) = q opt (2) in Table 1.

With constant τ but all other elements in q left free to vary, the optimization loop always stops around J(q)/J(q ref ) = 0.5. Loop termination is triggered by the | β(s)| ≤ 1 constraint being violated at zero frequency. Fig. 5(a) displays that within 7 iterations, the steepest descent algorithm with λ = 10 can lead to a reflected wave packet that has half the energy of that in the reference case. Smaller values of λ marginally decrease the final J, and the reflected p ′ (q opt , t) found with λ = 1 (leading to J(q opt )/J(q ref ) = .473) is shown in Fig. 5(b) along with p ′ ref and p ′ tar . The impedance and reflection properties of the TDIBC based on q ref and q opt are shown in Fig. 6, along with the data from Jones et al. [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF] that has been digitized from Fig. 11 in Monteghetti [START_REF] Monteghetti | Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized euler equations[END_REF]. Clearly, the optimization 2: Sensitivities of J computed at q = q ref through finite differences (FD) and automatic differentiation. Figure 5: Optimization at constant τ . (a) Convergence of J(q) towards a minimum value during the gradient descent optimization. For all three values of λ, the optimization algorithm is stopped due to the condition | β(s)| ≤ 1 being violated at zero frequency after {7, 71, 708} iterations for λ = {10, 1, 0.1}, respectively. (b) Reflected pressure perturbations: reference (q = q ref ), optimized (found with λ = 1, q = qopt from Table 1) and target p ′ tar .

process has modified the reactance properties of the TDIBC. Yet looking carefully at the differences between q ref and q opt in Table 1, it appears that the optimization at constant τ has essentially increased β ∞ and left the other elements of q ref unchanged in comparison. Partly because of this, partly because dJ/dβ ∞ is several orders of magnitude larger than the other sensitivities in Table 2, we explore how the optimization proceeds when we apply the further contraint of keeping β ∞ constant -on top of keeping τ constant. The result of the optimized q for constant τ and β ∞ is q opt2 , and it is displayed in Table 1. Clearly, the elements of q other than β ∞ and τ have now significantly changed. J(q opt2 )/J(q ref ) is decreased down to 0.55% as opposed to J(q opt )/J(q ref ) = 47.3%, proving that major improvements during the optimization can be gained by constraining β ∞ to remain constant. The decrease in J as the steepest descent algorithm proceeds 11 can be seen in Fig. 7(a), and the optimized reflected perturbation pressure is shown in Fig. 7(b). The impedance and reflection properties based on q opt2 are shown in Fig. 6. The optimization process used so far is carried out for an input wave packet with an amplitude spectrum centered around 1000 Hz. The optimized parameters in the TDIBC in Table 1 should correspond to an improved tubular ceramic liner that is tuned for the perturbation spectrum considered in Fig. 2. Two other broadband input spectra are considered next for the pressure perturbations, with f = 1500 Hz and f = 2000 Hz -see the caption of Fig. 2 for the exact dependence of p ′ on f . The input and the reflected perturbations at the probe location of the 1D impedance tube are shown in Fig. 8 for the three frequencies considered. The optimization procedure is applied to each input frequency, yielding the convergence properties displayed in Fig. 9. It can be observed that the optimization is able to proceed further for the two lower frequencies than for the highest one. The frequency-dependent optimized TDIBC settings, q opt3 , are shown in Table 3, and the corresponding impedance and reflection properties are gathered in Fig. 10. Finally, the original and optimized reflected perturbations are displayed in Fig. 11. It should be noted that in the case of the frequency-dependent optimization Figure 7: Optimization at constant τ and β∞. (a) Convergence of J(q) towards a minimum value during the gradient descent optimization. For λ = 10 8 , J(q)/J(q ref ) reaches 0.55% before the descent algorithm starts oscillating. We take q opt2 to be q recovered after 128 iterations. (b) Reflected pressure perturbations (reference and optimized), and target perturbation level p ′ tar .

leading to q opt3 , the target pressure perturbation p ′ tar is chosen to be the temporal mean of p ′ ref between t 1 and t 2 , as opposed to a constant value of p ′ ref (t 1 ). That is the reason why the q opt3 values on the left-hand column of Table 3 do not match those in the righthand column of Table 1, and why p ′ tar coincides with p ′ ref at the beginning in Figs. 5(b) and 7(b), but not in Fig. 11. One can see from Fig. 11 that p tar is somewhere in between p ′ ref at the beginning and end of the displayed time window. It is perhaps the appropriate point to illustrate that a pressure drift takes place, whereby the pressure in the tube decreases after the perturbations have gone past the probe location. The drift is most apparent at the highest input perturbation frequency, as can be seen in Fig. 11. This minor pressure drift is a well-known consequence of using a characteristic boundary condition which only balances out the perturbation field. A typical example is that of a non-reflecting boundary condition (B = 0) which will perfectly cancel-out the entering pressure wave but will not enforce any mean pressure value. That is, the computational domain receives no pressure feedback from the outside and the mean pressure field is free to depart from its initial state. is typically addressed by adding a linear relaxation term as in Poinsot and Lele [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF], at the expense of causing mild acoustic reflections as the boundary stiffness increases. Adding a linear relaxation term would modify the boundary's acoustic impedance profile, the accuracy of which is of paramount importance to the present study. Therefore, a minor pressure driftamounting to one-millionth of the bulk pressure -is deemed acceptable.

2D GIT flow

The setup described in subsection 3.2 is considered next. At the probe location (x = 70 cm and y = 5 cm in Fig. 4), the pressure perturbations are measured and used for the optimization procedure that seeks to bring them closest to a target level p ′ tar . The initial parameters determining the TDIBC behavior are those in q ref of Table 1. With these initial values, the plane wave packet can be seen in Fig. 12 coming from the inlet, while the solid portion corresponds to the perturbations reflected by the TDIBC operating with parameters set at q ref -see Table 1. the inlet towards the outlet as it is modified by the TDIBC. The pressure perturbations measured at the probe location with q ref are shown in Fig. 14(b) in red, where the amplitude of the perturbations has noticeably decreased from the initial 1 Pa value at parameter q opt3 (1000Hz) q opt3 (1500Hz) q opt3 (2000Hz) β∞ 5.000000000000000e-01 5.000000000000000e-01 the inlet. We note, in passing, the effect of the pressure drift discussed at the end of subsection 4.1, which is also visible in Fig. 14(b). The attenuation effect of the acoustic liner is most visible in Fig. 15, where the decrease in the sound pressure levels (SPL) measured along the wall opposite to the TDIBC is plotted. The initial TDIBC parameters in q ref correspond to the ceramic tubular liner (CT57) used in Jones et al. [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF], and leads to a SPL decrease of 9 dB under the present conditions: an input plane wave packet with frequency centered around f = 1500 Hz -see Fig. 2. The choice of this input frequency is justified in Section 5. The optimization procedure followed is identical to that described in the previous section, with constant τ and β ∞ . The cost function in Eq. ( 6) is computed between t 1 = 1.5 ms and t 2 = 10 ms, and p ′ tar is taken as constant and equal to the mean of p ′ ref during that time interval -see Fig. 14(b). The minimization of the cost function proceeds as shown in Fig. 14(a). After less than 30 iterations, the cost function goes below 10% of its original value. The spatio-temporal evolution of the plane wave packet obtained by using the optimized TDIBC parameters q opt,GIT of 4 is displayed in Fig. 13. The perturbations at the probe location and the sound pressure levels across the duct with the optimized TDIBC can be seen in Figs. 14(b) and 15, respectively. Finally, the reactance properties of the optimized liner are displayed in Fig. 16.

Figure 12: Plane wave packet evolution in the GIT test case, travelling at the speed of sound c = 344.37 m/s towards the right. The TDIBC is located at the bottom wall, between 20.3 cm ≤ x ≤ 60.9 cm -see Fig. 4. Simulation results obtained with the TDIBC parameters before optimization -q ref shown in Table 1.

Figure 13: Same as Fig. 12, but simulation results obtained with optimized parameters q opt,GIT with the reactance properties of Fig. 16.

Discussion

Before embarking on the discussion, it is convenient to note that a liner with fully absorptive power for perturbations of a single frequency f exhibits an impedance z such that Re (z) = 1 and Im (z) = 0 at frequency f . That concept proves useful when interpreting the changes induced by the optimization process on the impedance properties Figure 14: (a) Convergence of J(q) towards a minimum value during the gradient descent optimization. For λ = 100, J(q)/J(q ref ) reaches 6.9% while for λ = 1000, J(q)/J(q ref ) reaches 8%. In both cases, of the acoustic liner.

The results for the 1D impedance tube are discussed first. It is clear that the optimization process is able to proceed further after setting the bulk reflectivity β ∞ to a constant, yet it is interesting to observe how the acoustic properties of the TDIBC change when going from q ref to q opt in Table 1 by having acted almost exclusively on β ∞ . Fig. 6(b) shows that the frequency of the reflection peaks and troughs has not been altered, which is consistent with the poles of q opt remaining almost unchanged from those of q ref . However, the trough around 1 kHz has become shallower with the increase in β ∞ . The benefits of increasing β ∞ are less visible on the impedance data in Fig. 6, but a zoom into the curves in Fig. 6(c) close to 1 kHz (not shown) reveals that q opt brings the resistance closer to unity (from below) than q ref does, while the frequency of the zero parameter q ref q opt,GIT β∞ 5.000000000000000e-01 crossing towards 1 kHz in the reactance curves in Fig. 6(d) remains virtually unchanged.

Setting β ∞ to a constant allows the optimization to act on the oscillatory poles and their corresponding weights, leading to q opt2 in Table 1. The resulting parameters drastically improve the absorption properties of the TDIBC. Fig. 6(b) shows that the reflection trough around 1 kHz can be brought down to almost zero, consistent with a resistance closest to unity around that frequency in Fig. 6(c). The range of frequencies over which the reactance remains close to zero has been widened significantly with q opt2 , as observed from Fig. 6(d). Comparing p ′ ref with p ′ opt2 in Fig. 7(b) gives perhaps the most eloquent demonstration of the level of optimization that can be achieved when leaving τ and β ∞ constant for this input frequency centered around 1000 Hz. Keeping the optimization strategy unchanged, different input frequencies are analyzed next.

Fig. 9 shows that similar levels of optimization can be achieved with input wave packets that have amplitude spectra centred around 1000 Hz as well as 1500 Hz, yet for the 2000 Hz case the ability to optimize is severely hindered in comparison. The reason for this -at first surprising -behavior can be found in the constant value of τ . In the lossless limit, a tubular liner can be modelled as a quarter-wavelength resonator with resonant frequencies given by [START_REF] Monteghetti | Analysis and discretization of time-domain impedance boundary conditions in aeroacoustics[END_REF]

f n = (2n + 1)c 4l ( 8
)
where l is the cavity length. Since l and τ are related through l = cτ /2, the first two resonant frequencies obtained with the fixed τ = 4.79939 × 10 -4 s should be around 1042 Hz and 3124 Hz, which is consistent with the first two minima observed in the | β| figures based on q ref . The antiresonant frequency for this liner is located halfway in between the two absorption sweet points, i.e. close to 2000 Hz. It could therefore be anticipated that without changing the cavity length -and hence τ -, a drop in the optimization capacity should be observed when the frequencies to be attenuated are centred around 2000 Hz compared to 1000 Hz or 1500 Hz. A possible solution to this problem consists in shifting f n towards higher values by lowering l through τ . The 2D GIT flow is tackled in better conditions after the lessons learnt from the 1D impedance tube: 1) optimization at constant τ and β ∞ is most effective, and 2) frequencies close to 2000 Hz have limited optimization potential at the chosen τ . The CT57 liner is known to have an effectiveness at cancelling plane monochromatic waves which varies significantly with frequency. For a plane monochromatic wave of frequency 1000 Hz, Jones et al. [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF] report SPL drops of about 60 dB in the GIT experiment. At 1500 Hz, their data show a much deteriorated SPL drop of about 10 dB. An attempt is made here to find an alternative set of TDIBC parameters that improves upon the performance of the original liner, for frequencies where it is reasonable to expect that its effectiveness can be improved. For all these reasons, an input plane wave packet with amplitude spectrum centered around 1500 Hz is chosen: sufficiently far from the 1000 Hz optimum yet below the 2000 Hz limitations. From Fig. 14(a), it can be concluded that the optimization is able to proceed reasonably far, down to 6.9% of the initial value of the cost function. Although values below 1% are not reached as in the 1D impedance tube -see Fig. 9, it has to be emphasized that here all of the pressure perturbations do not come in contact with the TDIBC, which is only present on one of the two walls. Nevertheless, Fig. 15 illustrates that the original liner decreases the outlet sound pressure level by 9 dB compared to the inlet, and the optimized liner decreases them a further 12 dB. The impedance and reflection properties in Fig. 16 display once again the effect of the optimization. The reflection minimum in Fig. 16(b) has clearly been lowered and brought closer to the 1500 Hz range of frequencies.

The location of the probe at which the perturbations are minimized will necessarily have an effect on the resulting set of optimized TDIBC parameters. But the freedom to choose the probe location also allows one to be able to target more selectively where the perturbations ought to be minimized. Another possibility would be, of course, to embed into the cost function a pressure -or sound intensity -that is integrated throughout the outlet. It is only a matter of choice which does not imply any difficulty from an implementation point of view, since automatic differentiation can cope with any cost function definition as long as its expression can be differentiated.

The results shown so far are sufficient to illustrate that the optimization process can lead to different sets of optimized TDIBC parameters, even when starting from identical initial parameters q ref and keeping the input signal spectra unchanged. This can be seen by comparing, for instance, q opt3 at 1500 Hz in Table 3 with q opt,GIT in Table 4 -which is also obtained with the same input amplitude spectrum. Different flow conditions lead to different TDIBC optima, that is not surprising. What one may want to further assess, however, is how different two optima will be when all that is changed is the initial set of TDIBC parameters. That is assessed in Fig. 17 17. It appears that changing the initial set of values of q ref by up to 10% of their original value marginally modifies the optimum state, while the same test but with changes of up to 20% leads to an entirely different optimum. The conclusion drawn is that the initial q ref has a significant impact on the minimum that is found at the end of the optimization procedure, so that a poor choice of the initial q ref cannot, in general, be compensated by the optimization approach brought forward here. In the present study, the initial set of conditions is guided by the CT57 liner, which is already a fairly effective design that places the optimization procedure not too far away from optimum. In the work by Özkaya et al. [START_REF] Özkaya | Efficient design optimization of acoustic liners for engine noise reduction[END_REF], a first-stage, low-fidelity optimization is applied before the gradient descent algorithm, which is similar to placing the multiparameter space reasonably close to optimal values before diving in the search for an optimum point.

Conclusions

Numerical simulations are carried out using JAGUAR, a flow solver that features the possibility to apply a TDIBC [START_REF] Fiévet | Strong compact formalism for characteristic boundary conditions with discontinuous spectral methods[END_REF]. The parameters that control the behavior of the TDIBC are based on the oscillo-diffusive representation (ODR) of the acoustic reflection operator [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF][START_REF] Monteghetti | Energy analysis and discretization of nonlinear impedance boundary conditions for the time-domain linearized euler equations[END_REF], with values chosen to reproduce accurately the properties of the ceramic tubular liner CT57 [START_REF] Jones | Benchmark data for evaluation of aeroacoustic propagation codes with grazing flow[END_REF]. The TDIBC parameters are modified in order to optimize the behavior of the TDIBC at dampening pressure perturbations with different broadband frequency contents. The optimization follows a gradient descent algorithm, where the gradients of the cost function with respect to the ODR parameters has been obtained by means of the automatic differentiation compatibility of JAGUAR [START_REF] Cardesa | Adjoint computations by algorithmic differentiation of a parallel solver for time-dependent pdes[END_REF].

A 1D impedance tube test case is considered first, which reveals that optimizing the TDIBC at constant bulk reflectivity β ∞ and time delay τ proves most effective -as 21 opposed to letting these two parameters vary. Radically improved dampening properties compared to the original TDIBC are obtained following this strategy. The effectiveness of the optimization is found to be limited when the input perturbations come close to the anti-resonant frequency of the tubular liner modelled by the TDIBC. An obvious workaround to this limitation is to carry out the optimization by keeping τ constant yet fixed to a smaller value. Finally, a 2D test case of the GIT is chosen to study a more realistic flow configuration, where the liner is known to be operating at off-design conditions. The optimization succeeds at adapting the set of TDIBC parameters to recover an improved absorptive performance. Throughout the present work, the optimization is applied to quantities appearing in the ODR of the acoustic liner's reflection operator. For liners characterized by few and simple geometric quantities, there are mathematical models that can approximate the corresponding ODR coefficients describing the reflection or impedance properties. The converse is not as trivial, since not all impedance properties can be realized using a geometry-based mathematical model. As a consequence, it is not guaranteed that the optimized parameters found in the present work can be generated with a modified CT liner. Since the bulk reflectivity and the tube length are fixed, the remaining free parameter that can change the CT liner impedance is the tube diameter alone. That will not suffice to select the values of the 12 parameters left to adjust within the parameter vector q. Additional geometric parameters would be needed in the mathematical model, which would for instance involve tubes of more complex cross-section. As a result, improvements to the present approach include using the geometry-based mathematical model inside the flow solver, with routines that translate the geometrical parameters directly into ODR coefficients. In this manner, the optimization could be carried out based on updates of the input geometric parameters rather than on the ODR coefficients, thus limiting the optimization space to the realm of readily obtainable acoustic liners.

Figure 1 :

 1 Figure 1: Sketch of the 1D impedance tube test case. The inlet packet of pressure fluctuations travels towards positive x at the sound speed c (top). At x = L = 1.8 m, a TDIBC is located which interacts with the incident wave packet, reflecting it back towards the inlet (bottom).

Figure 2 :

 2 Figure 2: (a) Time series of the pressure perturbation of the incident wave packet signal and (a) its corresponding single-sided amplitude spectrum. The input wave packet for a chosen frequency f is given by p′ = A sin (2πf (t -3/f )) exp -σ 2 f (t -3/f )2 . For all cases considered in the present study,

Figure 3 :

 3 Figure 3: Pressure perturbations p ′ measured about the ambient value pa = 101325 Pa at the probe location x = 0.3 m.The oscillations for t < t 1 are due to the incident wave packet characterized in Fig.2travelling from the inlet at x = 0 towards the TDIBC wall at x = 1.8 m. The interval between t 1 and t 2 is chosen to study the properties of the pressure perturbations reflected by the TDIBC back towards the inlet.

Figure 4 :

 4 Figure 4: Sketch of the grazing incidence tube (GIT) test case. L = 812 mm, H = 50 mm. The right running pressure waves are damped under the effect of the TDIBC on the lower wall.

Figure 6 :

 6 Figure 6: (a) Argument and (b) modulus of the reflection coefficient β(s) based on q ref , qopt and q opt2in Table1. (c) Real and (d) imaginary part of the impedance z(s) based on q ref , qopt and q opt2 . z = (1 + β)/(1 -β).

Figure 8 :

 8 Figure 8: Pressure perturbations p ′ measured about the ambient value pa = 101325 Pa at the probe location x = 0.3 m in the 1D impedance tube.Dashed lines correspond to the right-running wave packet coming from the inlet, while the solid portion corresponds to the perturbations reflected by the TDIBC operating with parameters set at q ref -see Table1.

Figure 9 :

 9 Figure9: Convergence of J(q) towards a minimum value during the gradient descent optimization at three different input wave packet frequencies -and hence three different J(q ref ). In all cases, λ = 100. For f = 2000 Hz, the optimization ended soon due to the | β(s)| ≤ 1 constraint being violated at 0 Hz. This is not the case for the other two frequencies, where computations are terminated manually after vanishing returns upon further iterations.

Figure 10 :

 10 Figure 10: (a) Argument and (b) modulus of the reflection coefficient β(s) based on q ref and q opt3 for the three analysed frequencies f . (c) Real and (d) imaginary part of the impedance z(s). z = (1+ β)/(1-β).

Figure 11 :

 11 Figure 11: Reflected pressure perturbations (reference and optimized), and target pressure level p ′ tar . (a) f = 1000 Hz, (b) f = 1500 Hz, (c), f = Hz.

Figure 15 :

 15 Figure14:(a) Convergence of J(q) towards a minimum value during the gradient descent optimization. For λ = 100, J(q)/J(q ref ) reaches 6.9% while for λ = 1000, J(q)/J(q ref ) reaches 8%. In both cases, the optimization is terminated due to the | β(s)| ≤ 1 condition being violated at 0Hz. (b) Pressure fluctuations (reference and optimized) at the probe location (x = 70cm, y = 5cm) for the 2D GIT case, and target pressure level p ′ tar .

Figure 16 :

 16 Figure 16: (a) Argument and (b) modulus of the reflection coefficient β(s) based on q ref and q opt,GIT . (c) Real and (d) imaginary part of the impedance z(s). z = (1 + β)/(1 -β).

  . The different data have been obtained by changing each element of q ref -except for β ∞ and τ -through multiplication by (1 + θ), where θ is chosen from a uniform random distribution spanning the range -|q ref (m)|R ≤ θ ≤ |q ref (m)|R. For R = 0.1 and R = 0.2, the convergence towards a minimum is plotted in Fig.

Figure 17 :

 17 Figure17: Convergence of J(q) towards a minimum value during the gradient descent optimization, with f = 1000Hz and λ = 100 in all cases. Two convergence curves are obtained by perturbing each element of q ref by a random factor of up to ±0.1 (or ±0.2) of its value.

Table 1 :

 1 . In the GIT experiment, t 1 = 1.5 ms and t 2 = 10 ms. Input quantities parametrizing the TDIBC properties for the case where Nν = 0 and Ns = 4 in Eq. (

	While p ′ tar (t) can in general be a function of time, we choose it to remain fixed so that p ′ tar (t ≥ t 1 ) = p ′ ref (t = t 1 ). The choice is entirely driven by our goal to reduce

Table 3 :

 3 q opt3 found by optimizing at constant τ and β∞ for frequency-dependent p ′ tar .

	5.000000000000000e-01

Table 4 :

 4 q opt,GIT is found by optimizing at constant τ and β∞ in the 2D GIT case. Input plane wave packet centered around 1500Hz. Found with λ = 100 after 32 iterations.

	5.000000000000000e-01

European Center for Research and Advanced Training in Scientific Computation

For clarity, all complex numbers are written with a tilde, while exclusively real-valued ones are not.

Acknowledgments

This work was part of the 3C2T project, which was financially supported by the STAE foundation (Science and Technology for Aeronautics and Space) through grant number CDT-R050-L00-T00. The project involved the French aerospace laboratory (ONERA) through the internal project JEROBOAM, the European Center for Research and Advanced Training in Scientific Computation (CERFACS) and the Institute of Fluid Mechanics of Toulouse (IMFT) under the supervision of Christophe Airiau (IMFT). J.I.C. acknowledges funding from the People Program (Marie Curie Actions) of the European Union's Seventh Framework Program (FP7/2007-2013) under REA grant agreement n. PCOFUND-GA-2013-609102, through the PRESTIGE program coordinated by Campus France. J.I.C. is grateful to R. Roncen for enlightening discussions on interpreting impedance and reflection data.


Marcel Bouché

N Bottinelli

M Hedde

P Jouquet

Y Capowiez
email: yvan.capowiez@inra.fr

An explicit definition of earthworm ecological categories

Keywords: classification, ecological traits, functional group, endogeic, epigeic, anecic, epi-anecic

In the early 70s, the soil biologist Marcel Bouché classified French earthworms species into defined ecological categories. This classification system was immensely successful and is still widely used to describe earthworm functional groups even outside of Europe. Bouché used morpho-anatomical traits to differentiate three main categories: epigeic, anecic and endogeic. However, the way species are assigned to a category was not explicitly described in Bouché's work. Thus, some earthworm species can still be assigned to two categories depending on the way researchers interpret Bouché's description. To solve these issues and avoid unnecessary controversies, we applied PCA and random forest models to the seminal data of Marcel Bouché (earthworm morpho-anatomical traits). Their assignment to Bouché's three main categories allowed us to statistically redefine the different categories and determine which traits are the most influential. We found that the three main traits were skin pigmentation (from none to black), body length (mean of the minimal and maximal values) and skin coloration (yes or no), followed by 10 other morphological and anatomical traits. We then used this approach to assign a likely category to all of the species studied by Bouché, resulting in a new triangular graph including other categories such as epi-anecic, endo-anecic, epiendogeic and intermediate. Finally, we calculated the percentage that each species belongs to each main ecological category. This represents a paradigm shift and may change our vision of earthworm communities enabling the computation of the percentage of anecic, endogeic and epigeic species at the community level and thus overcoming the limits and debate about fixed ecological categories for each species.

Introduction

Marcel Bouché defined a widely used earthworm classification based on three ecological categories (epigeic, anecic and endogeic) at the beginning of the 1970s. The concept was first introduced in a chapter of an edited book [START_REF] Bouché | Relation entre les structures spatiales et fonctionnelles des écosystèmes illustrés par le rôle pédobiologique des vers de terre[END_REF], further explained in detail in his book [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] and finally summarized again in a conference paper a few years later [START_REF] Bouché | Soil organisms as components of ecosystems[END_REF]. This last paper was particularly successful since it described, for the first time, a triangle defined by three poles corresponding to three main ecological categories and positioned some typical earthworm species within this triangle. Despite the fact that these reports were all published in French, they had incredible success and the categories were widely adopted, at least for Lumbricidae species. As a consequence, the two main references [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF][START_REF] Bouché | Soil organisms as components of ecosystems[END_REF] have been highly cited and remain consistently cited today: in the last 50 years, about 10% of the papers on earthworms have cited at least one of these two references (Fig. S1).

Marcel Bouché was especially influenced, among others, by the seminal work of [START_REF] Wilcke | Über die vertikale Verteilung der Lumbriciden im Boden. Einteilung der Lumbricidenfauna nach Michaelsen[END_REF] and suggested that three main criteria can be used to classify earthworms, their: (i) vertical distribution, (ii) ecology and (iii) morphology and with obvious, although not quantified, relationships between these criteria.

However, because vertical distribution can vary depending on the environment, climate and thus the season, Marcel Bouché only developed his classification using the two last criteria and first defined ten morpho-physiological (1972) and later fifteen morpho-anatomical (1977) characteristics (Table S1). The second criteria (i.e. the ecology) remained, however, unquantified and relied on Bouché's personal observations or experience (1972). He also clearly stated that some earthworm species could be classified between the poles of the triangular graph. Thus four additional classes were defined and termed 'epi-endogeic ', 'epi-anecic', 'endo-anecic' and 'intermediate' (the latter class means that the species is positioned between the three main classes). Examples of such species are Microscolex dubius (Fletcher, 1887) for epi-endogeic, Lumbricus terrestris (Linnaeus, 1758) for epi-anecic, Aporrectodea caliginosa meridionalis [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] for endo-anecic and Allolobophora chlorotica (Savigny, 1826) as intermediate.

So, why revisit this concept? Firstly, because some species studied by Marcel Bouché in his book were not assigned to an ecological category (something that Marcel Bouché did not comment upon). Secondly, because despite the impressive success of this classification and the countless services this classification has rendered to scientists trying to understanding soil ecology, some recurrent controversies have plagued earthworm literature, e.g., the status of Aporrectodea trapezoides Duges, 1828[START_REF] Baldivieso-Freitas | Earthworm abundance response to conservation agriculture practices in organic arable farming under Mediterranean climate[END_REF] or Allolobophora chlorotica [START_REF] Couteulx | Burrowing and casting activities of three endogeic earthworm species affected by organic matter location[END_REF]. These inconsistencies are difficult to prevent since no general rules were given regarding the use of the 10 to 15 characteristics of [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF][START_REF] Bouché | Soil organisms as components of ecosystems[END_REF] and their relative importance. For species that were classified as occupying an intermediate position between two or three poles, it is even more difficult since, in this case, Bouché did not provide a formal definition. Futhermore, although the triangle (1977) is well known to soil ecologists, the method used to generate it is now ignored by most of us (how many morpho-ecological characteristics were used?; were they all useful?) and has even been forgotten by Marcel Bouché himself (Bouché, personal communication).

Hence, to revisit the ecological categories, we first carried out a literature review looking for the most studied earthworm species worldwide to determine whether their ecological category was unambiguously defined and largely accepted in scientific publications. Then we used both the traits described in Bouché's book (1972) and the ecological category assigned to the species, to study the most influential characteristics and the minimal set of characteristics that could be used to assign an earthworm species to an ecological category. Finally, all the species mentioned by Bouché were projected onto a definitive triangular graph. This also means that, for each species, our approach goes beyond a fixed ecological category (this use was not encouraged by Bouché) and instead computes the percentages by which they belong to the three main ecological categories. Our hypothesis was that this numerical approach, based on morpho-anatomical data, would provide a more accurate classification of earthworm ecological categories.

Material and Methods

Bibliography review

To define the most currently studied earthworm species at the global scale, we benefited from the recent 11 th International Symposium of Earthworm Ecology, held in Shanghai in the summer of 2018. Using the book of Abstracts (http://www.isee11.org/Home/Menu/30), we recorded every earthworm species cited in oral and poster presentations. We only retained species that were cited at least twice and obtained an initial list of 36 species. Since our main focus was soil ecology, Eisenia fetida and E. andrei were removed from the list, since these species are mainly found in composts and mainly studied in ecotoxicology and vermicompost production. For each species in the initial list, we then searched the Web of Science database (http://apps.webofknowledge.com) publications using the earthworm genus and species names as keywords. We applied the following rules to select the manuscripts: 5 years old (2013)(2014)(2015)(2016)(2017)(2018) and if less than 100 manuscripts were found for that period, we considered older manuscripts until 100 manuscripts were selected (when possible). We then retained species for which at least 10 manuscripts could be found, obtaining a final list of 26 species. For each selected publication, when whole access was possible, we searched 'endogeic', 'epigeic' and 'anecic' as keywords, so that the three intermediate categories were also detected.

Thus we determined whether the species under consideration was assigned to an ecological category (only the ecological category 'intermediate' was not detected).

When available, we also noted the explanations given by the authors to explain this assignment.

Definition of ecological categories based on earthworm morpho-anatomical traits

Forty morpho-anatomical traits measured on 125 species were recorded in Bouché's book (1972) even if Bouché used only 10 and then later 15 of them to define the ecological categories. This dataset is available from the Dryad Digital Repository (http://doi.org/10.5061/dryad.g7046). Not all the species defined by Bouché are currently accepted by the scientific community (see for example http://drilobase.org/), but we decided to use all the data available in the book in order to replicate how Bouché defined his ecological categories. As a first important rule, sexual traits were not taken into account in our study because [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] never mentioned them in the definition of the ecological categories and we feared that this would have given too much weight to phylogeny. Our dataset (Table S2) included 125 species (18 epigeic, 25 anecic, 20 endogeic, 6 epi-anecic, 10 epiendogeic, 2 endo-anecic, 5 intermediate and 39 unspecified) characterized by 23 traits (11 qualitative and 12 quantitative) potentially relevant for distinguishing ecological categories. To build a model capable of assigning an earthworm species to an ecological category, we adopted the following procedure: first, we selected the morpho-anatomical traits which were the most important drivers of the differences among the three ecological categories (epigeic, anecic and endogeic) considering that these are the three main poles from which the other four categories were derived. To do this, we carried out a random forest classification to identify the most important predictors characterizing the three groups among the 23 traits. The importance of each predictor was determined by assessing the increase in the mean square error (MSE) between observations and out-of-bag predictions. The significance of the importance of each predictor on the three groups was tested with a permutation method involving 500 replicates of each tree, using the 'rfPermute' package in R. The selected traits were then used to build a PCA model using the 'PCAmixdata' package in R. From this PCA model, we added other species, which were not assigned to a category but described by [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] as supplementary individuals (6 epi-anecic, 10 epi-endogeic, 2 endo-anecic, 5 intermediate and 39 unspecified). The triangle was defined by the three barycenters of each main ecological category in the two first axes of the PCA (the barycenters were computed only using species assigned to each category by [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF]. Once the three poles of the triangle were defined, we computed the coordinates of each species as the distance to each pole in the two first axis of the PCA. These distances were first used to determine whether a species was within or outside of the triangle. If within, its coordinates were simply given by the distances to the three poles. If outside, then the species was projected onto the closest pole or side of the triangle and its coordinates determined by recomputing the distances to the three poles after the projection. The triangle was further divided into seven sub-regions to arbitrarily define the seven ecological categories. We set the limits of the main categories to 75% to get true archetypes and then cut the central space into four regions by setting the limits to 60%.

Results

Assignments to an ecological category found in the literature review

Among the 36 species first selected, we found more than 10 available publications for only 26 of them. More than 50% of these species were assigned to only one ecological category (Table 1). Seven species were assigned to two ecological categories. Among these, A. trapezoides and Metaphire guillelmi were assigned to two of the main categories (anecic-endogeic-epigeic). The remaining five were assigned to an intermediate category and this was the case for the well-known L. terrestris which is mainly cited as an anecic but sometimes classified as an epianecic (Table 1). Next, two species were classified in three categories: A. longa was mainly cited as anecic (84%) but also as endogeic (6%) and endo-anecic (10%) whereas L. rubellus was mainly classified as epigeic (72%) but also epi-anecic (2%) and epi-endogeic (25%). Lastly, two species were classified in four categories: Amynthas corticis was stated either endogeic (13%), epigeic (37%), epi-anecic (13%) or epi-endogeic (37%), whereas Amynthas gracilis was considered as anecic (14%), endogenic (14%), epigeic (43%) or epi-endogeic (29%). In general, in the above-cited studies the authors did not provide data or other evidence to support the assignment of the species to an ecological category. Indeed, the percentage of articles with such supporting information varied from 0 to 43% depending on the species under consideration (Table 1). It is worth noting that authors usually referred to earthworm behavior such as feeding (37.2% of explanations), followed by their vertical distribution in the soil profile (25.5%) and gallery organization (28.3%). In contrast, data related to earthworm morphology and anatomy were rarely used since earthworm pigmentation and size represented only 2.8 and 6.2% of the data, respectively (data not shown).

Characteristics and minimal dataset for the definition of ecological categories

The random forest classification model performed very well with regards to small out-of-bag errors (4.8%). Among the 23 traits used, only 13 were considered important for differentiating the three main ecological categories (Table S4). The most important traits included skin coloration (qualitative trait), body length (quantitative trait) and skin pigmentation (qualitative trait) (Table S3). This emphasized the importance of skin attributes in the definition of ecological categories. The rest of the traits which discriminated the categories were, in order of importance: body weight, the antero-posterior gradient, the dorso-ventral gradient, diameter, epithelium rigidity, muscle structure, typhlosole type, flattening index, transversal and longitudinal furrows (definitions can be found in Table S4). Other evaluated traits had a negligible effect on the differences found among categories.

The PCA model generated with the 13 important traits and 63 species (those assigned to one of the three main categories by [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] showed a good separation of the epigeic, anecic and endogeic species (Fig. 1) although intra ecological group variability was rather large. The first axis (22% of total variance explained) clearly separated the three ecological categories whereas the second axis (13% of total variance explained) opposed epigeic species to anecic and endogeic species. The projection of supplementary species (those unassigned or those assigned to ecological categories other than the three main ones) on the two axes of the PCA model showed that the 6 epi-anecic species were located between epigeic and anecic species, whereas the 10 epi-endogeic species overlapped with endogeic species. The five intermediate and the two endo-anecic species had variable positions in the PCA. Several species which [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] did not previously assign, preferentially overlapped with the endogeic species.

Definition of ecological categories

The scores for each species obtained from the PCA were projected on a ternary plot (Fig. 2) and some well-studied species and species mentioned on Bouché's triangle (1977) are highlighted. We observed a relatively good agreement between Bouché's triangle (1977) and our triangle. For each species, the percentage belonging to the three main ecological categories and their assignment to one of the seven ecological categories are presented in Table S4. Our model indicated that among the 125 species, endogeic species dominated (31%), followed by epigeic (17%), epi-anecic (18%) and anecic (14%). The other categories were scantly represented with less than 10% for each (Table 2). The classification accuracy between our model and Bouché's groups (1972) showed high values for epi-anecic (100%), epigeic (82%) and endogeic (70%) medium values for intermediate (60%) and anecic (52%) and very low values for epi-endogeic (25%) and endo-anecic (0%). The 13 traits characterizing the four most abundant ecological categories used in this study (epigeic, anecic, endogeic and epi-anecic) are presented in Table 3. From this analysis we determined a set of characteristics for each category.

Epigeic species are small sized and light in weight. They are all pigmented, usually red with a gradient along the body. The caudal part of the body is sometimes flattened and the epithelium is soft. The muscle structure is pinnate. They do not have a longitudinal furrow and few have transversal furrows. The shape of the typhlosole is mostly massive. Anecic species are larger and heavier than epigeics.

They are all pigmented, with a brown color and a gradient along the body. The clitellum and the caudal part of the body can be flattened and the epithelium is rigid.

The muscle structure is mostly pinnate. They do have longitudinal and transversal furrows. The shape of the typhlosole is only pinnate. Endogeic species are of medium size and weight. They are all non-pigmented, thus their coloration (visible through their skin) is mostly pink with no gradient along the body. The clitellum can be flattened and the epithelium is soft. The muscle structure is mostly elementary and radial. They all have longitudinal furrows but few have transversal furrows. The shape of the typhlosole is mostly bilamelated. Epi-anecics share similar traits with anecics except the absence of longitudinal furrows, which is a trait shared with the epigeic species.

Discussion

Accurate and unambiguous definitions of the earthworm ecological categories

Marcel Bouché's ecological categories have been outstandingly successful and there is no doubt that they have been instrumental in unifying studies on earthworm ecology. However, our literature review showed that, even for some very common Lumbricidae such as A. chlorotica, L. terrestris, and A. nocturna, it can be challenging to clearly assign ecological categories It should be noted however that confusion regarding the classifications were generally between a main category (anecic, endogeic, epigeic) and one of the sub-categories. This reflects the fact that for some authors [START_REF] Hoeffner | Feeding behaviour of epi-anecic earthworm species and their impacts on soil microbial communities[END_REF], these ecological categories are indeed hierarchical: the earthworm species is first anecic and then further described as epianecic, for example (as often observed for L. terrestris; van Groenigen et al. 2019). This is in contrast with Bouché's opinion that these categories were rather exclusive as showed by his definition of the 'intermediate' ecological category (in between anecic, endogeic and epigeic categories).

There are also difficulties with earthworm species from Asia, with the striking examples of two species (A. gracilis and A. corticis) assigned to four different ecological categories. We assume that the lack of clear rules impeded the authors from assigning these species to one ecological category. Moreover, the ecological categories of traits defined mainly on Lumbricidae in Europe may not perhaps be fully applicable to earthworms from other families and continents. This makes it clear to us that the time has come to revisit these ecological categories, not to destroy them but to rebuild their foundations in order to apply the approach to either (i) difficult cases or (ii) other continents (if applicable). This is also supported by the availability of modern statistical methods that were out of reach at the beginning of the 1970s (random forest approaches).

Our method successfully reconstructed and accurately defined the three main ecological categories and four additional categories. Every new species for which the 13 traits are measured or estimated can thus be positioned in the triangle and assigned an ecological category (assuming that this classification made mainly on Lumbricidae can be extended to other families). Our approach could also be used to confirm the use of these categories for non-lumbricidae earthworms from other parts of the world such as Asia. The hierarchization of the factors (given by the random forest classification) to define the ecological categories was not unexpected.

However, the three main characteristics are very basic traits (skin coloration, body length and skin pigmentation) and can be easily assessed on adult earthworms. The remaining nine characteristics do include some internal anatomy (typhlosole, muscle structure) requiring dissection or the use of categorical characteristics (flattening index and epithelium rigidity) requiring expertise or the establishment of clear descriptions and boundaries. Regarding the last criteria (epithelium rigidity), it is worth mentioning that Briones and Alvarez-Oterao (2018) proposed earthworm body wall thickness as a method for inferring membership to an ecological category.

Despite the general agreement between our classification and Bouché's seminal work there were some discrepancies. The most striking differences concern the anecic vs epi-anecic categories. In brief, in [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF], anecics combined the brown-headed large earthworms with the giant Mediterranean earthworms (mainly of the Scherotheca genus), and the epi-anecic earthworms are mainly red-headed large earthworms of the Lumbricus genus. In contrast, in our study anecics only included most of the giant Mediterranean earthworms whereas the red and brownheaded large earthworms from the Lumbricus and Nicodrilus genera were epi-anecic (with a few exceptions). The main difference between these two categories in our classification was mean size and the existence of transversal furrows (Table 3). This may be seen as non-optimal or illogical for some researchers, especially those that define epi-anecics using behavioral traits [START_REF] Hoeffner | Feeding behaviour of epi-anecic earthworm species and their impacts on soil microbial communities[END_REF]) something Marcel Bouché never did. This suggests that, using multivariate analysis, the more extreme anecic earthworms (giant Mediterranean worms) represent the norm for anecic earthworms. Although this observation may hold true for statistical reasons only, in section 4.2, we explain how this may not be a problem. Besides these discrepancies, we also observed clear agreements with Bouché's work. The most impressive is the position of A. chlorotica which we found in the middle of the triangle and was thus assigned to the 'intermediate' category (in the middle of the three main categories), exactly how Bouché predicted in his book (1972, p. 269): « A. chlorotica is morphologically assigned to an intermediate position between the three categories, even if the green morph has manners more epigeic than the leucotypic morph ». This is striking since no papers were found that cited or acknowledged this strange position for this current species. In our literature review, A. chlorotica was always assigned to the endogeic category with the exception of only one report (Le [START_REF] Couteulx | Burrowing and casting activities of three endogeic earthworm species affected by organic matter location[END_REF] where it was assigned to the epi-endogeic category. In our opinion, this agreement reinforces our classification and the seminal work and intuition of Marcel Bouché who never explained in his book why this species was given an intermediate position: it was so obvious to him that no further explanations was needed. However, our study makes these statements fully explicit and also fully available to non-French speakers.

From fixed categories to the use of fuzzy logic

The use of sole categories to classify earthworm species was originally criticized by [START_REF] Bouché | Soil organisms as components of ecosystems[END_REF] himself who advocated that anecic, epigeic and endogeic are the three main poles of a continuous distribution of earthworm species between these three main ecological strategies. However, despite its success in terms of scientific citations, Bouché's seminal work has been frequently misapplied and most, if not all, the researchers involved in earthworm ecology, including the authors of the present manuscript, have used the ecological categories as categorical data. It is time for a rediscovery of what Bouché had already stated in the 1970s: earthworm species representing a typical earthworm category are rare. In our study, we found that only 16 species out of the 125 belonged 100% to only one ecological category. These were: two Dendrobaena and one Dendrodrilus as epigeic ; most of the Prosellodrilus, Hemigatsrodrilus monicae, Ailoscolex lacteospumosus and Allolobophora graffi and A. zicsii as endogeic ; and Nicodrilus gognus as anecic. It is once again striking that only one species was 100% anecic, however most of the Scherotheca species were more than 90%.

Slicing of the triangle based on arbitrary thresholds may have induced bias and could be a problem for assigning species close to these borders. However, following the idea of Marcel Bouché, we believe it is time for a paradigm shift regarding these categories and we should now accept that each earthworm species is indeed defined by three percentages, one for each main ecological category (see Table S4). This is not a drastic increase in complexity of the definition of these categories and could help solve recurrent controversies regarding the assignment of an earthworm species to a unique and fixed ecological category [START_REF] Bastardie | A radio-labelled study of earthworm behaviour in artificial soil cores in term of ecological types[END_REF], Felten and Emerling 2009[START_REF] Baldivieso-Freitas | Earthworm abundance response to conservation agriculture practices in organic arable farming under Mediterranean climate[END_REF], especially the species not assigned by Marcel Bouché. In this study, we have defined, using automatic borders inside the triangle, the four other categories that are used but which were not accurately defined by [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF]. We also provide some information regarding their main characteristics (excluding the two less common categories) to provide support for further discussion on this topic. We do not believe that our proposed separation is definitive. It is worth noting that all the misclassified species using our methodology were part of the mixed categories (epianecic, endo-anecic, epi-endogeic and intermediate; Table 2). In contrast, we advocate for the wider use of percentage of belonging to the three main (well accepted and defined) categories.

The perspectives of our work go beyond dispute resolution, it opens the avenue for new kinds of approaches in earthworm ecology. Traditionally, earthworm data are first analyzed using total abundance and biomass and then using ecological categories and their putative reaction to different stresses such as tillage, pesticide application (and thus exposure) or soil compaction. We believe our approach can be used almost unambiguously to compute the percentage of the anecic, endogeic and epigeic categories in whole earthworm communities based on fuzzy logic estimates of their membership to the ecological categories (i.e. the three percentages). This approach is analogous to the Community Weighted Means widely applied in trait-based ecology [START_REF] Hedde | Functional traits of soil invertebrates as indicators for exposure to soil disturbance[END_REF].

We also believe that this is a first important step in the study of relationships between earthworm communities and important soil functions associated with their activities [START_REF] Frazao | Soil structure formation and organic matter distribution as affected by earthworm species interaction and crop residues placement[END_REF][START_REF] Van Groeningen | How fertile are earthworm casts? A meta-analysis[END_REF]. However, this goal will only be achieved if we develop new traits that can be defined for earthworm species, particularly traits linked to the ecology and behavior of species and not only those based on morpho-anatomical characteristics.

Conclusions and prospects for the future

We re-analyzed the data of [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] to redefine the three main ecological categories for earthworms. By doing this, we were also able to compute the percentage by which each species belongs to each category based on data of 13 morpho-anatomical traits. We assume this will change the way soil ecologists will analyze earthworm communities in the future.

One of our aims was also to propose a shared procedure for computing this percentage for species worldwide. It is not possible to provide a universal spreadsheet for doing this since in the PCA we used mixed categorical and continuous data, which does not give access to the computations of the projections on the two main axis. Instead we now propose that if research groups send us their data (13 traits) collected for a species not already described in Bouché's book (1972), we will analyze them and provide the percentages of belonging to the three main categories. We will also share the results with the wider scientist community (upon agreement with the authors) on websites such as 'BETSI' or 'Drilobase'. Since the traits could be measured or estimated in different ways by different researchers, our percentages will be computed on several datasets (for example Pontoscolex coeruthrus (Muller, 1856) from different parts of the world) enabling the computation of robust mean values.

Figure and Table captions

Fig. 1: Projection of the 125 earthworm species, including non-assigned species and those belonging to four additional ecological categories, described by [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] in the PCA defined using the 13 more important morpho-anatomical traits (assessed by a random forest approach using only the three main ecological categories).

Fig.2:

Ternary plot defined by the three main categories (anecic, endogeic, epigeic) where the 125 species were projected. Seven zones were defined to assign each species to one ecological category (but see Table S4), some typical species of each category are listed.

Table 1: Summary of the literature review for the 27 species selected: (i) the number of articles studied for each species; (ii) the assignment of the species in ecological categories; (iii) the percentage of articles where the assignment to an ecological category was justified and (iv) the number of articles where ecological categories were not mentioned. [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] with the model-based approach used in the present study.

Table 3: Summary of the morpho-anatomical traits characterizing the most frequent ecological categories studied by [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF].

(ii) the assignment of the species in ecological categories; (iii) the percentage of articles where the assignment to an ecological category was justified and (iv) the number of articles where ecological categories were not mentioned. [START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF] with the model-based approach used in the present study.

Epigeic

Table 2 :

 2 A comparison of the ecological categories assigned by

Table 2 :

 2 number of studies in which the eartworm species was not assigned to an ecological category by the authors A comparison of the ecological categories assigned by

	Species	Total	Anecic	Endogeic	Epigeic	Epi-anecic	Epi-endogeic	Endo-anecic	Intermediate	Justification	No mention*
						(%)					(number)
	One category
	Aporrectodea caliginosa	86		100						28	19
	Aporrectodea icterica	19		100						37	2
	Octolasion lacteum	19		100						5	7
	Octolasion cyaneum	33		100						6	9
	Aporrectodea rosea	64		100						11	29
	Octolasion tyrtaeum	36		100						8	9
	Pontoscolex corethrurus	52		100						8	9
	Amynthas robustus			100
	Lumbricus castaneus	25			100					16	7
	Eisenia veneta	52			100					6	10
	Dendrobaena octaedra	52			100					2	4
						1

*

Table 3 :

 3 Summary of the morpho-anatomical traits characterizing the most frequent ecological categories studied by[START_REF] Bouché | Lombriciens de France. Ecologie et Systématique[END_REF].

	Morpho-anatomical traits	Epigeic	Endogeic	Anecic	Epi-anecic
	Length (mm)	62	79	349	143
	Diameter (mm)	2	3	9	6
	Weight (mg)	518	751	13235	3054
	Pigmentation	100%	0%	100%	100%
	Skin coloration	Red (86%)	Pink (61%)	Grey (100%)	Grey (55%)
	Antero-posterior gradient	100%	0%	100%	100%
	Dorso-ventral gradient	95%	0%	84%	100%
	Flattening index	No (43%), Tail (29%)	No (45%), Clittelum (29%)	Clitellum-Tail (47%), Tail (47%) Clitellum-Tail (45%), Tail (23%)
	Epithelium rigidity	38%	0%	89%	91%
	Muscle structure	Pinnate (95%)	Elementary (47%), Radial (16%)	Pinnate (63%)	Pinnate (95%)
	Longitudinal furrow	0%	34%	89%	5%
	Transversal furrow	43%	100%	100%	95%
	Typhlosole type	Massive (57%)	Bilamelated (58%)	Pinnate (100%)	Pinnate (68%)

Acknowledgments

We thank colleagues from the TEBIS network (http://www.reseau-tebis.fr/) for fruitful discussions about ecological traits and their use.