Jérémy Sadet

Franck Massa

Thierry Tison

Isabelle Turpin

Bertrand Lallemand

El-Ghazali Talbi

Homotopy perturbation technique for improving solutions of large quadratic eigenvalue problems: Application to friction-induced vibration

Keywords: Quadratic Eigenvalue Problem, Homotopy, Perturbation, Projection, Nonlinear vibration, Friction-induced vibration

This paper puts forward a projection technique for accurately calculating solutions of large Quadratic Eigenvalue Problem. The aim here is to stabilize the complex eigensolutions whilst reducing residual errors, especially when considering significant damping contribution or asymmetric stiffness matrices. Hence, more confident results can be obtained in the frequency band of interest. To achieve this, high order modes, calculated using the homotopy perturbation technique, are introduced in the projection step of the classical method. This numerical proposal is a generalization of the classical projection, based only on normal modes of the associated undamped problem. To evaluate the efficiency of the suggested method, a finite element application dedicated to a friction-induced vibration problem is investigated.

Introduction

In mechanical engineering, the modal analysis of large undamped structures discretized by the Finite Element Method (FEM) is performed using the Generalized Eigenvalue Problem (GEP) from the definition of mass and stiffness matrices. However, when damping is of main interest or in other applications, such as vibroacoustics, fluid dynamics or friction-induced vibration problems, solutions are the results of a Quadratic Eigenvalue Problem (QEP) for which a complete and practical survey can be found in [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. In this case, large discrete models require a projection onto a relevant subspace. In modal analysis, the classical method uses a basis built with a set of eigenvectors of the associated undamped model. For the complex applications cited above, this basis is not generally optimal and can lead to significant errors. For example, for a frictioninduced vibration problem, one of the reasons for potential errors is that the coupled normal modes are a poor representation of the model dynamics subjected to friction. Many methods (sub structuring, component mode synthesis, condensation, reanalysis techniques, etc.) are based on the construction of reduced models by means of projection on a basis of displacements. It is known that these methods provide accurate results when properly used, or in other words, when the reduced subspace spanned by the columns of the rectangular projection matrix represents a close approximation of the full order model [START_REF] Balmès | Parametric families of reduced finite element models. theory and applications[END_REF]. It is also known that increasing the subspace size with additional vectors of the same type slowly improves the solutions and is not interesting for numerical cost reasons.

It is thus more interesting to investigate other categories of subspaces already used in other contexts such as Taylor series [START_REF] Massa | Structural modal reanalysis methods using homotopy perturbation and projection techniques[END_REF], perturbations techniques [START_REF] Chen | Comparison of several eigenvalue reanalysis methods for modified structures[END_REF][START_REF] Kirsch | Nonlinear dynamic reanalysis of structures by combined approximations[END_REF][START_REF] Ma | Combined approximation for reanalysis of complex eigenvalues[END_REF][START_REF] Jian-Jun | Structural modal reanalysis for large, simultaneous and multiple type modifications[END_REF], Padé approximant [START_REF] Damil | An iterative method based upon Padé approximants[END_REF][START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF] or Homotopy Perturbation Method (HPM) introduced by He [START_REF] He | Homotopy perturbation technique[END_REF]. At first, the HPM was used in the literature to calculate nonlinear solutions for different applications. Duigou et al. [START_REF] Duigou | Iterative algorithms for non-linear eigenvalue problems. application to vibrations of viscoelastic shells[END_REF] have associated homotopy, asymptotic numerical techniques and Padé approximants to investigate the vibrations of damped sandwich structures whereas Boumediene et al. [START_REF] Boumediene | A reduction model for eigensolutions of damped viscoelastic sandwich structures[END_REF] have developed a reduction method for solving the complex nonlinear eigenvalue problem regarding viscoelastic structures. Sun et al. [START_REF] Sun | A homotopy method for determining the eigenvalues of locally or non-locally reacting acoustic liners in flow ducts[END_REF] have provided numerical solutions for soundpropagation computations using the homotopy method. Next, Lee et al. [START_REF] Lee | Natural frequencies of thin rectangular plates using homotopy-perturbation method[END_REF] have searched for asymptotic solutions of nonlinear problems with the HPM to yield natural frequencies. For vibroacoustic applications, Claude et al. [START_REF] Claude | Eigensolutions to a vibroacoustic interior coupled problem with a perturbation method[END_REF] put forward a numerical method to solve non-symmetric eigenvalue problems and discussed of the efficiency of eigensolvers depending on Taylor or asymptotic expansions in [START_REF] Claude | Numerical comparison of eigenvalue algorithms for vibroacoustic problems[END_REF]. Recently, the same authors have suggested computing damped eigenfrequencies and modes with a high order Newton solver, based on homotopy and perturbation techniques [START_REF] Claude | Study of damped vibrations of a vibroacoustic interior problem with viscoelastic sandwich structure using a high order newton solver[END_REF]. Secondly, several authors have considered the HPM to estimate modified solutions of a perturbed problem by reanalysis or by Reduced Order Model (ROM). Indeed, Sliva et al. [START_REF] Sliva | A study of the eigenvalue sensitivity by homotopy and perturbation methods[END_REF] suggested solving the new eigenvalue problem with the combination of a homotopy transformation and the perturbation method. Next, Li et al. [START_REF] Li | Modal modification of damped asymmetric systems without using the left eigenvectors[END_REF] have worked on the calculation of complex eigenpairs of modified asymmetric systems as a way of creating a Campbell diagram of the modified rotor system. Massa et al. [START_REF] Massa | Multi-level homotopy perturbation and projection techniques for the reanalysis of quadratic eigenvalue problems: The application of stability analysis[END_REF] have integrated two-level homotopy techniques to estimate instabilities of a friction-induced vibration problem. Do et al. [START_REF] Do | A global strategy for the stability analysis of friction induced vibration problem with parameter variations[END_REF] used the HPM to reanalyse modal bases of modified components of a mechanical system before using them in the component modal synthesis. More recently, Massa et al. [START_REF] Massa | From homotopy perturbation technique to reduced order model for multiparametric modal analysis of large finite element models[END_REF] have defined a ROM for a linear modal analysis by considering a projection basis made of high order perturbed eigenvectors.

In this paper, we suggest using the HPM to compute QEP high order modes. These modes are used to build an enriched projection basis before solving the reduced QEP resulting from a dynamic analysis of structures discretized by the FEM. With the classical method, variations of the complex eigensolutions values are often observed when the number of normal modes in the projection basis is modified. With the suggested method, the aim is to stabilize the results of the QEP while reducing the residual errors. The quality of the eigensolutions is thus improved and results are proven to be more accurate.

Section 2 reviews the QEP and its classical way of resolution. Section 3 generalizes the classical projection basis to higher order eigenvectors. Next, section 4 presents an example of a friction-induced vibration problem. Finally, section 5 offers a summary of our conclusions for this novel strategy.

Quadratic Eigenvalue Problem definition

Although some problems in structural dynamics are defined by non-symmetric matrices, only the right quadratic eigenvalue problem will be discussed in the sequel. Let us consider the QEP given by Eq.(1) :

K + s i C + s 2 i M ψ i = 0 (1)
where s i and ψ i , are respectively the i th eigenvalue and the i th right eigenvector.

M, C and K are respectively the [N × N] mass, damping and stiffness matrices of the discretized mechanical system, with N representing the number of Degree Of Freedom (DOF). The case where the three matrices are Hermitian positive definite, or positive semidefinite for C and K, corresponds to a structure modeled with viscous damping. This QEP corresponds also to a damped mechanical system where components are coupled by normal contact entries. Eigenvalues of such problems are real or complex conjugate and come in pairs. For gyroscopic systems, M and K are Hermitian, with M definite positive, and C skew Hermitian (C = -C H , with C H the conjugate transpose of C). In this case, eigenvalues are purely imaginary. When the mechanical system is subjected to friction (such as friction-induced vibration problems), tangential loads induced a non-symmetric stiffness matrix resulting from the tangential coupling. Even if damping is not considered, complex solutions of the latter QEP can reveal positive real parts which make it possible to determine the stability of the system. Other types of problems such as vibroacoustic, fluid mechanics studies, control theory in the field of automatic and even numerical optimization can also be reduced to a QEP. For a complete description of these applications, the interested readers can refer to the paper of Tisseur et al. [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. The suggested algorithm can be applied to the aforementioned problems. In this paper, we will specifically illustrate the case of a friction-induced vibration problem because the friction effect can introduce large errors in the evaluation of complex eigensolutions. The well-known global QZ algorithm [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF] is efficient for solving the QEP of small discrete numerical systems (typically when N<1000). For large Finite Element (FE) models, a projection onto a T subspace is needed to build a reduced approximate modal problem that is compatible with the QZ algorithm. In practice, the subspace T is made up of the first normal modes of the associated conservative system given by Eq.( 2) :

K -ω 2 i M φ i = 0 (2)
where ω i and φ i , are respectively the i th natural circular frequency and the i th normal mode.

Considering both the projection basis and the linearization of the QEP (according to the definition in [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]), one can obtain a GEP of the second companion form given by Eq.(3).

Au i = s i Bu i (3)
where

A = -T T KT 0 0 T T MT , B = T T CT T T MT T T MT 0 , u i = q i s * i q i
and ψ * i = Tq i , s * i and ψ * i are approximations of the QEP eigensolutions in Eq.(1).

The above strategy is common and largely implemented in commercial finite element softwares. In most of the practical problems, particularly the modal analysis with viscous damping, solutions of Eq.( 3) are very satisfactory. When studying the stability of rubbing systems, some discrepancies arise which could alter the judgment of the behaviour of the mechanical system. These errors can be reduced by adding appropriate modes calculated from the HPM presented in the next section.

Homotopy resolution of the QEP equation

Many physical phenomena are modelled using systems of nonlinear differential equations, which is a straightforward way to describe their behaviours. As these equations are generally difficult to solve, many alternative and powerful methods have been developed over the last few years, for instance, the decomposition method [START_REF] Cherruault | Decomposition methods: A new proof of convergence[END_REF], the variational iteration method [START_REF] He | Variational iteration method -a kind of non-linear analytical technique: some examples[END_REF][START_REF] Wang | Variational iteration method for solving integro-differential equations[END_REF], the HPM [START_REF] He | Homotopy perturbation technique[END_REF], Generalized Homotopy Method (GHM) [START_REF] Vazquez-Leal | Generalized homotopy method for solving nonlinear differential equations[END_REF] and the differential transform method [START_REF] Chen | Solving partial differential equations by two-dimensional differential transform method[END_REF]. Each method is limited to a special class of integro-differential equations and short reviews are available for different applications [START_REF] He | A short review on analytical methods for a fully fourth-order nonlinear integral boundary value problem with fractal derivatives[END_REF][START_REF] He | A short review on analytical methods for the capillary oscillator in a nanoscale deformable tube[END_REF]. The Homotopy Perturbation Method, used in this paper to build the projection space, has been improved from 2000 to present days [START_REF] He | A coupling method of a homotopy technique and a perturbation technique for non-linear problems[END_REF][START_REF] He | New interpretation of homotopy perturbation method[END_REF][START_REF] He | Homotopy perturbation method with two expanding parameters[END_REF] and applied to study different problems integrating nonlinear oscillators [START_REF] Yu | Homotopy perturbation method with an auxiliary parameter for nonlinear oscillators[END_REF][START_REF] He | The simpler, the better: Analytical methods for nonlinear oscillators and fractional oscillators[END_REF]. The following section details an implementation of the homotopy perturbation technique for solving a QEP problem where matrices can be asymmetric.

Homotopy Perturbation technique

Solutions of non-linear equations can be obtained using the HPM which combines the standard homotopy in topology and the perturbation technique. The idea behind the HPM is to separate the linear and nonlinear parts of the problem and to introduce an embedding parameter ε to highlight the non-linearity. Assuming that the embedding parameter is small, approximate solutions of the non-linear equations can be written as a power series in ε.

Considering the general case of the QEP, the basic idea of the suggested methodology is to introduce the HPM in two ways. The first one concerns the symmetric and asymmetric contributions of the mass and stiffness matrices. The decomposition in Eq.(4) represents the symmetric contribution, the S index, and the asymmetric contribution, the A index. Decomposition is presented here for both mass and stiffness matrices for the sake of generality but the asymmetric contribution can easily be suppressed in the case of a purely symmetric matrix.

M = M S + εM A K = K S + εK A (4)
where M S , M A , K S , K A stand for the symmetric and asymmetric decomposition of the mass and stiffness matrices, respectively.

The second decomposition is related to damping which is completely considered as a non-linear part of Eq.( 1). The decomposition of the damping matrix C is then different as all entries are considered to be asymmetric contributions such as :

C = εC A (5)
According to the HPM, the i th eigensolution of Eq.( 1) is developed as a power series in ε with a nominal term, the (0) index, and higher order contributions, indexed by (n), such as :

s i = s (0) i + εs (1) i + . . . + ε n s (n) i + . . . + ε d s (d) i ψ i = ψ i (0) + εψ (1) i + . . . + ε n ψ (n) i + . . . + ε d ψ (d) i ( 6
)
where d is truncation order.

Substituting Eq.( 4), Eq.( 5) and Eq.( 6) into Eq.( 1), the QEP is rewritten as follows :

(K S + εK A ) + s (0) i + εs (1) i + . . . + ε n s (n) i + . . . + ε d s (d) i εC A + s (0) i + εs (1) i + . . . + ε n s (n) i + . . . + ε d s (d) i 2 (M S + εM A ) ψ (0) i + εψ (1) i + . . . + ε n ψ (n) i + . . . + ε d ψ (d) i = 0 (7)
From Eq.( 7), it is possible to identify the n th homotopy order and the high order vectors associated with that order. The next subsections describe the form of the nominal (0 th order), the 1 st order and any n th order of the high order vectors.

Identification of the nominal terms

Terms of the 0 th order are simply obtained setting ε = 0 in Eq.( 7) leading to a GEP of the form :

K S + s (0) 2 i M S ψ (0) i = 0 (8)
where s

(0) i is the i th nominal eigenvalue and ψ (0) i is the i th nominal eigenvector.

Eq.( 8) has no asymmetric contributions. Moreover, substituting s (with j, the imaginary unit), the problem is equivalent to the linear eigenvalue problem in Eq.( 2) and can be computed with a classical eigenvalue algorithm, for example a Lanczos solver, with eigenvectors normalized according to the mass matrix such as :

ψ (0)T i M S ψ (0) i = 1 (9)

Identification of the 1 st order terms

Considering only the 1 st order terms, Eq.( 7) is rewritten as followed :

K S + s (0) 2 i M S ψ (1)
i + 2s

(0) i M S ψ (0) i s (1) i + G (1) i = 0 G (1) i = K A + s (0) i C A + s (0) 2 i M A ψ (0) i (10)
Premultiplying ( 10) by ψ (0) T i and considering Eq.( 8), the 1 st order eigenvalue is written as :

s (1) i = -ψ (0) T i G (1) i 2s (0) i (11)
Eq.( 10) can be rewritten to be of the form of Eq.( 8) with a non-zero right member :

K S + s (0) i M S ψ (1) i = F (1) i ( 12
)
where

F (1) i = -2s (0) i M S ψ (0) i s (1) i -G (1) i Unfortunately, K S + s (0) 2 i M S is not invertible. Another equation is needed to determine ψ (1)
i .

The chosen solution follows Wang's proposal for the calculation of eigenvector derivatives in [START_REF] Wang | Improved approximate methods for computing eigenvector derivatives in structural dynamics[END_REF]. The idea is to express each perturbed eigenvector as a linear combination of the 0 th order eigenvectors ψ (0) i , i.e normal modes, of the system. For the 1 st order eigenvector, the linear combination reads:

ψ (1) i = N k=1 c (1) k ψ (0) k (13)
In practice, the number m of available normal modes is lower than N. A common solution is to complete the development with a static correction term Γ

i . The linear combination Eq.( 13) is then rewritten as:

ψ (1) i = m k=1 c (1) k ψ (0) k + Γ (1) i Γ (1) i = K -1 S F (1)
i -

m k=1 ψ (0) k F (1) i w (0) 2 k ψ (0) i ( 14
)
For i = k, coefficients c (1)
k are calculated with Eq.( 15).

c

(1)

k =i = ψ (0) T i F (1) i w (0) 2 k -w (0) 2 i ( 15
)
If i = k, Eq.( 15) cannot be evaluated. The coefficients c

(1)

k=i have to be determined by other means. A second convenient normalization, given by Eq.( 16), is considered. Instead of directly considering a perturbed mass normalization homogeneous to Eq.( 9), this equation links the 0 th order eigenvector ψ (0) i and the high order eigenvector ψ i to the full mass matrix. The aim is to limit the complexity of the following developments and use the orthonormal property of Eq. [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF].

ψ (0) T i Mψ i = 1 (16)
Introducing Eq.( 4) and Eq.( 6) in Eq.( 16) and keeping only the 1 st order terms, Eq.( 17) reads as :

ψ (0) T i M S ψ (1)
i + ψ

(0) T i M A ψ (0) i = 0 (17)
Substituting the 1 st order eigenvector by the linear combination given by Eq.( 13) and considering that normal modes are M S orthonormals, the c

(1) k=i coefficient can be calculated with Eq. [START_REF] Sliva | A study of the eigenvalue sensitivity by homotopy and perturbation methods[END_REF].

c (1) k=i = -ψ (0) T i M A ψ (0) i (18)

Identification of the n th order terms

Following the developments detailed in section 3.3, expressions to evaluate any n th order of the i th eigensolutions can be determined and are given by Eq.( 19) to Eq. [START_REF] He | Variational iteration method -a kind of non-linear analytical technique: some examples[END_REF].

K S + s (0) 2 i M S ψ (n) i = F (n) i (19)
where

F (n) i = -2s (0) i M S ψ (0) i s (n) i -G (n) i (20) G (n) i = K A ψ (n-1) i + C A n-1 k=0 ψ (k) i s (n-k-1) i + M A n-1 k=0 n-1-k l=0 s (l) i s (n-k-l-1) i ψ (k) i + 2M S   n-1 k=1 s (0) s (k) i ψ (n-k) i + 1 {n≥3} n 2 -mod((n+1),2) l=1 n-l k=l+1 s (l) i s (k) i ψ (n-k-l) i   + M S   n 2 k=1 s (k) 2 i ψ (n-2k) i   (21)
where mod(x, 2) is the modulo 2 of x, x corresponds to the greatest integer less than or equal to x, also denoted floor(x) and where the indicator function 1 {n≥p} equals 1 if n is superior to p and 0 otherwise.

s (n) i = -ψ (0) T i G (n) i 2s (0) i (22) ψ (n) i = m k=1 c (n) k ψ (0) k + Γ (n) i Γ (n) i = K -1 S F (n) i - m k=1 ψ (0) T k F (n) i ω (0) 2 k ψ (0) i (23) c (n) k =i = ψ (0) T i F (j) i ω (0) 2 k -ω (0) 2 i (24) c (n) k=i = -ψ (0) T i M A ψ (n-1) i (25)
The computation of perturbed eigensolutions has to solve a linear system of size N (Eq.( 14) or Eq.( 23)) for each mode which can be time consuming for large systems. However the symmetric contribution of the stiffness matrix involved in these expressions can be decomposed once to minimize the additional computational cost.

Construction of the projection basis

With the classical method, the projection basis T is built with nm normal modes of Eq.( 2), such as :

T = [φ 1 . . . φ nm ] (26)
Since this basis is obtained with the associated conservative system's solutions, the subspace contains no information of damping nor any matrix asymmetric contributions. Thus, the key idea of the suggested method is to integrate damping and asymmetric information into the projection basis with additional order modes calculated from HPM. In a sense, the suggested method, further named HOPEP in this paper, retains the same characteristics of the classical method for large problems, namely linearization of the QEP, projection to reduce the corresponding GEP and then resolution with a global algorithm. HOPEP can then be considered as a generalization of the classical method.

Consequently, the projection basis T in HOPEP takes the following form :

T = φ (0) 1 . . . φ (0) nm ψ (1) 1 . . . ψ (d) m ( 27
)
where m is the number of complex perturbed modes at each order. It is worth noting that the same number of normal modes is retained in expansions Eq.( 14) and Eq.( 23). This solution guarantees a constant number of normal modes for the computation of the perturbed modes during the iterations of the HOPEP algorithm. Several tests have also shown that increasing the number of normal modes did not significantly improve the results.

When dealing with complex values in a numerical analysis, it is common practice to separate the real and imaginary parts of mathematical expressions. This solution allows the preservation of as much information as possible but doubles the number of modes in the projection basis, and therefore the time to solve the QEP. A more economic solution considers only the real parts of perturbed eigenvectors. It is important however to keep sufficient information in the real part before decomposing complex quantities. It is a common practice [START_REF] Fuellekrug | Computation of real normal modes from complex eigenvectors[END_REF] to scale the complex matrix with the maximum value in each column and retain only the modified real part.

As the projection basis T in HOPEP is composed of a set of modes from different types or orders, it is mandatory to proceed with orthonormalization phase. Here, two numerical methods have been identified. The first method is the so-called iterative Gram-Schmidt process [START_REF] Jamai | Influence of iterated Gram-Schmidt orthonormalization in the asymptotic numerical method[END_REF]. This method is very accurate and allows for the computation of an orthormal basis preserving the initial size. The second method, implemented in the Structural Dynamics Toolbox [START_REF] Balmès | Structural Dynamics Toolbox (for use with MATLAB)[END_REF], orthonormalizes the projection basis using a preconditioned Cholesky decomposition and eliminates possible collinear vectors using a singular value decomposition. An advantage of this second solution is a lower numerical cost. For both methods, vectors of the basis are normalized thanks to the identity matrix.

Once the different calculations detailed above have been performed, the modified projection basis T is used in Eq.( 3) to solve approximate solutions of the QEP.

Many tests have shown that HOPEP's performance relies on a good balance between the number of normal modes, the number of perturbed eigenvectors and the maximum homotopy order. The flowchart in Fig. 1 depicts a strategy to efficiently build the projection basis. The algorithm mainly focuses on the sequence and criteria used to define the number of normal modes and the maximum homotopy order. Indeed, the number of perturbed eigenvectors (m) is fixed to the number of mode of interest since a stabilization of all complex eigensolutions is required in the entire frequency band. It is worth noting that every change in the projection basis is followed by a new evaluation of the QEP. Starting with several normal modes equals to the number of modes in the frequency band of interest, additional normal modes or homotopy orders are iteratively added to the projection basis. The process is achieved by monitoring the evolution of the residual error and the convergence of eigensolutions. For each mode, the residual error is calculated using the Manhattan norm of the residual vector defined by Eq.( 28):

R i = K + s * i C + s * 2 i M ψ * i 1 (28)
It is expected that perturbed eigenvectors of any order added to the projection basis, will reduce the residual errors for all complex modes compared to those calculated at order 0. However, in practice, it is more relevant to first increase the number of normal modes rather than homotopy order if the residual error is not reduced for all modes, or at least a large number of them. It is important to note that an increase in normal modes has no effect on perturbed eigenvectors previously calculated since they required only the first m normal modes in Eq.( 14) or Eq. [START_REF] Moler | An algorithm for generalized matrix eigenvalue problems[END_REF]. Perturbed eigenvectors already calculated can thus be stored and reused. The convergence criterion is evaluated using the relative errors between eigensolutions of two successive orders. In practice, it is sufficient to monitor only the real parts of eigensolutions. The algorithm stops when relative errors for all modes of interest are inferior to a given threshold, otherwise a new homotopy order is added to the projection basis.

In the next section, the different projection basis configurations detailed above are studied in a finite element application dedicated to a friction-induced vibration problem.

Application of the HOPEP method

The stability of a brake system can be analyzed in the Lyapunov sense through the frequency method. This method involves performing a Complex Eigenvalue Analysis (CEA), in which the system of equations is linearised around a stationary state. Positive real parts of complex eigenvalues allow to state for unstable modes.

All tests performed on the brake model described in section 4.1 require a starting phase presented in section 4.2. Sections 4.3 and 4.4 analyse the effects of key parameters thanks to the convergence of complex eigensolutions and to the evolution of the residual error. The algorithm presented in section 3.5 is applied in section 4.5. In each case, the frequency band of interest is inferior to 20 kHz. In this frequency band, 74 complex modes are computed from the classical projection method.

For all tests, homotopy developments are limited to a 6 th order which represents a good balance between precision and numerical cost. Results are analysed thanks to the evolution of the QEP residues and relative errors on complex eigensolutions. For complex eigenvalues, standard formulations compare the real parts Eq.( 29) and the natural circular frequencies Eq.( 30), separately. The chosen relative error in complex eigenvectors, defined by Eq. [START_REF] He | A short review on analytical methods for the capillary oscillator in a nanoscale deformable tube[END_REF], is of the most accurate since major errors in only one vector component are noticed. It is normal to consider that a value below 10% denotes a very good agreement. An eigenvector correlation is performed with the most widely used Modal Assurance Criterion (MAC), as recalled in Eq.( 32) for complex quantities.

e Re(si) = Re(s * (n+1) i ) -Re(s * (n) i ) Re(s * (n+1) i ) (29)
e |si| = |s * (n+1) i | -|s * (n) i | |s * (n+1) i | ( 30
)
e ψ i = ψ * (n+1) i -ψ * (n) i ψ * (n+1) i + ψ * (n) i (31)
M AC ik = ψ * (n+1)H i ψ * (n) k 2 ψ * (n+1)H i ψ * (n+1) i ψ * (n)H k ψ * (n) k ( 32
)
Figure 2: FE model of the simplified brake system.

Studied model

The case study is a simplified brake system composed of two pads and a disc. An exploded view of the brake FE model and a detailed view of a pad are depicted in Fig. 2.

Each pad is constituted of a lining with a centered diagonal slot and a back plate. The full mesh is composed of 23302 nodes and 19104 linear brick elements. The coupled model contains 66490 DOF. Material properties are cast iron for the disc, steel for the backplates and a sintered material for the linings.

Material behaviour of the disc and backplates is supposed isotropic but transverse isotropic material is considered for the linings. The dimensions and material properties of each component are reported in Table 1.

A Rayleigh damping has been chosen to model the damping of the three components. Mass and stiffness damping coefficients, namely α and β, as well as friction coefficient µ between pads and disc are reported in Table 2.

To complete the model description, fixed boundary conditions for all directions are applied on interior diameter nodes of the disc (the red markers in Fig. 2) and all translations along x and y are fixed for the pads. Moreover, a pressure of 25 bars, represented by the red arrows, is applied on the external pad surfaces and a rotation velocity of 8.17 rad/s around the z axis is applied to the disc to initiate sliding friction.

Preliminary computations

The FE model of the brake system is pre-computed in the commercial software Abaqus. The stiffness K 0 , damping C 0 and mass M 0 matrices of the uncoupled system (i.e bloc diagonal matrices) are exported in a dedicated software developed in Matlab. To initiate HOPEP calculations, symmetric and asymmetric decompositions of the matrices are computed using normal and tangential connectivity matrices T N and T f respectively together with Eq. [START_REF] He | New interpretation of homotopy perturbation method[END_REF] to Eq.( 37) :

M S = T T N M 0 T N ( 33
)
M A = T T f (µ)M 0 T N -M S ( 34
)
K S = T T N K 0 T N ( 35
)
K A = T T f (µ)K 0 T N -K S ( 36
)
C A = T T f (µ)C 0 T N (37)

Influence of the projection basis arrangement

This first test shows the mutual influence of the number of normal modes and the homotopy order for the convergence of complex eigensolutions. In order to do this, three configurations of the projection basis are considered:

• a gradual increase of sets of normal modes (from nm to 7×nm)

• nm normal modes completed by sets of m perturbed eigenvectors up to order 6

• 2×nm normal modes completed by sets of m perturbed eigenvectors up to order 6

For the three configurations, values of parameters nm and m are fixed to the number of complex modes of interest (nm=m=74). Perturbed eigenvectors are decomposed into real and imaginary parts and the projection basis is orthonormalized using the iterative Gram-Schmidt process in order to preserve all available information. Fig. 3 shows the convergence evolutions of the eigensolutions namely, real part of eigenvalues, natural circular frequencies and eigenvector norms, as a function of two successive orders. Residual evolution is represented as a function of orders. It is worth noting that, for the first projection basis configuration, an order stands for a set of normal modes. As we would like to observe an improvement for all results in the frequency band of interest, convergence and residual evolutions are monitored by respectively the maximum relative errors and the maximum residual errors which are calculated for all modes.

The four graphs in Fig. 3 clearly show that the first configuration of the projection basis does not lead to improved results although more than five hundred normal modes are kept at the last order. The same remark applies to the second configuration. Perturbed eigenvectors have no effect when they are associated with the minimal number of normal modes. For the third projection basis configuration, where the number of normal modes is doubled, results are very interesting. In fact, in all cases the evolution of the three convergence criteria present the same decrease, attesting to the stabilization of solutions. The main improvement is observed for real parts and eigenvector norms where maximum relative errors calculated between order 0 and order 1 are relatively high denoting large discrepancies for some modes. Increasing the homotopy order significantly reduce these errors which can reach very low values for real parts. For example, reductions of two orders of magnitude are obtained with only two homotopy orders. For circular natural frequencies the maximum relative error is already low (below 1%) between order 0 and order 1 denoting a good evaluation of these quantities with the classical projection basis. The evolution of residual errors follows the same trend and confirms a gradual improvement in the QEP solving with increasing homotopy order. In conclusion, these first results give us two main characteristics of the suggested method. Firstly, the use of perturbed eigenvectors efficiently improves the QEP solutions provided that the size of the normal modes basis is greater than the number of mode of interest. Secondly, the stabilization of QEP solutions can be monitored with only the real part of eigenvalues.

Influence of normalization scheme

The second test aims to study the influences of the orthonormalization scheme as well as the decomposition of complex perturbed eigenvectors using either real and imaginary parts (labelled R+I) or modified real parts (labelled Rm). The projection basis can be orthonormalized either with the Gram-Schmidt process (labelled GS) or with the Cholesky decomposition (labelled Chol). In all four cases, the configuration of the projection bases follows the best out of the previous test, namely a normal mode basis of 148 eigenvectors completed with sets of 74 perturbed eigenvectors up to order 6.

As before, Fig. 4 shows the maximum relative errors of the eigensolutions as a function of two successive homotopy orders and the maximum residual errors as a function of homotopy orders. Regarding the stabilization of eigensolutions, there is no significant difference between the four configurations meaning that the convergence of complex eigensolutions is not affected by the orthonormalization and the decomposition schemes. Regarding the evolution of the residual errors, the global trend is the same for all configurations with a reduction of errors as a function of homotopy order. The best configuration remains the one which uses the Gram-Schmidt process on the real and imaginary parts of perturbed eigenvectors. The main difference deals with the gap in residual errors values mainly due to the choice of the modified real part decomposition and in a lesser extent, to the Cholesky decomposition. However, in both cases the gap compared to the best configuration does not exceed half an order of magnitude. It is thus interesting to consider the modified real part decomposition for lower numerical cost reason. Indeed, in that case, the size of the projection basis is half the size of the other decomposition scheme.

Application of the HOPEP algorithm

In the HOPEP algorithm, the iterative construction of the projection basis relies on the residue criterion and the stabilization criterion which govern an increase of respectively the number of normal modes and of the homotopy order. The residue criterion at order n is validated if the difference between the residual error at that order and the residual error at order 0 is negative. Ideally, this criterion should be validated for all complex modes. In practice, a reasonable threshold can be fixed to 95% of the number of modes. With the second criterion, solutions of the QEP are considered stabilized if relative errors between real parts of complex eigenvalues for two successive orders are inferior to a given threshold, typically 5%.

In this section, the behaviour of these two criteria is studied through three configurations of the projection basis. The first configuration corresponds to the initialization phase of the algorithm where the size of the normal mode basis is fixed to the nm=74 complex modes of interest. For the second configuration, the normal mode basis is chosen in a frequency band increased by 10% (up to 22 kHz) and contains nm=93 modes. In the same way, the frequency band for the third configuration is increased by 20% which provides a normal mode basis with nm=113 modes. For the three configurations the projection bases are progressively completed by the modified real part of perturbed eigenvectors up to order 6 and are orthonormalized thanks to the Gram-Schmidt process. We are interested by the 74 first complex modes of the QEP and the threshold for the stabilization criterion is fixed to 5%. The results of the first configuration are reported in Fig. 5 with the progression of the homotopy order (from top to bottom). A red dot indicates that the residual improvement or the stabilization is not valid for the related mode. With this configuration, the normal mode basis is not sufficient to stabilize the complex modes in the higher frequency band and the residual error is not improved for 10% of modes, no matter what the homotopy order is.

In the same way, the results of the second configuration are presented in Fig. 6 and in Fig. 7 for the third configuration. In both cases, validations of the two criteria have notably improved. However, Fig. 6(b) shows that six homotopy orders are required to stabilize modes in the entire frequency band. We can also note that residual errors are improved for the majority of modes at each step of the iterative process. As expected with the third configuration, results again improved faster since three homotopy orders are enough to stabilize all modes.

Table 3 presents a focus on real parts and frequencies calculated with the three configurations from order 0 (normal modes only) up to order 3 in a way that shows the behaviour of the results according to the setting of the projection basis more specifically. In this table, positive values of real parts, which correspond to unstable modes, are highlighted in light gray. Apart from mode number 46 where results are stable, we can globally observe once again that real parts of eigenvalues are mainly affected, both in signs and in values. Important variations of real part values can be observed for modes number 47, 51 and 52 but they stabilized quickly with homotopy orders in the third configuration (nm=113).

For mode number 53 the real part is always negative for the first configuration (nm=74) indicating a stable mode. For the second (nm=93) and the third configurations, this real part becomes positive as soon as one homotopy order is taken into account but remains negative if only normal modes are considered. We can also observe that real parts and frequencies completely stabilize with three homotopy orders for the third configuration.

With regard to modes 48 to 50, important fluctuations of real parts are observed depending on the number of normal modes used in the projection basis. Regarding the first configuration, starting with three stable modes, two unstable modes are finally detected if homotopy orders are added. Nevertheless, large variations of the values of the real parts lead to a doubtful perception of this result. On the other hand, for the second and the third configurations, the state of the modes is completely different since only the mode number 49 is detected as unstable. This situation clearly shows that the stability of the eigensolutions is governed by a good balance of normal modes and homotopy orders. Indeed, with the highest number of normal modes without homotopy order, modes number 48 and 50 would be considered as unstable. However, the amount of normal modes remains essential to improving the convergence of the stabilization as can be observed for the third configuration compared to the second one.

All observations given above are also strengthened by correlations and relative eigenvector errors reported in Table 4. Comparison between eigenvectors of the same number can be arguable because it can be seen that some MAC values, for modes 47 to 49, do not correspond to any agreement between modes shapes. There is no reference here and eigensolutions are then analyzed in the sense of convergence. One can observe a very good improvement in correlations and a drastic decrease in errors as soon as both normal mode number and homotopy orders increase. This is particularly true for the third configuration at order 3 where relative eigenvector errors are below 1% which denotes quasi identical vectors.

This test shows that the numerical cost can be reduced if the normal mode basis is increased first. This condition is controlled by the residual improvement at the first step of the algorithm and needs a reduction in residue for a large number of modes, 98% in that case. Moreover, if less normal modes are selected, simulations in the second configuration tell us that a better quality of complex eigensolutions can be obtained with higher homotopy orders. Consequently, for an application such as friction-induced vibration, the mode's stability status can be better evaluated.

Conclusion

In this study, an homotopy projection basis has been suggested as a generalization of the classical projection basis used to solve a QEP. The key idea with the HPM is here to consider the asymmetric parts of the mass and stiffness matrices and full damping matrix as non-linear parts of the problem. The development of perturbed eigensolutions is formulated in a way to reduce the numerical cost. An algorithm is also suggested to efficiently build a sufficient projection basis, again, minimizing the additional numerical cost. The numerical analyses clearly highlight that the new projection basis allows for the good stabilization of complex eigensolutions in the entire frequency band. Moreover, a significant reduction of QEP residue confirms the confidence of the results.

Figure 1 :

 1 Figure 1: Construction of the projection basis.

Figure 3 :

 3 Figure 3: Influence of projection basis content. Evolution of errors as a function of homotopy orders for: (a) real parts; (b) natural circular frequencies; (c) complex eigenvector; (d) QEP residual.

Figure 4 :

 4 Figure 4: Influence of projection basis orthonormalization. Evolution of errors as a function of homotopy orders for: (a) real parts; (b) natural circular frequencies; (c) complex eigenvectors; (d) QEP residual.

Figure 5 :

 5 Figure 5: (a) Residual criterion checking per mode and (b) stabilization criterion checking per mode as a function of homotopy orders for the first configuration (nm =74). A red dot indicates an unsatisfactory check.

Figure 6 :Figure 7 :

 67 Figure 6: (a) Residual criterion checking per mode and (b) stabilization criterion checking per mode as a function of homotopy orders for the second configuration (nm =93). A red dot indicates an unsatisfactory check.

Table 1

 1 Geometrical and material parameters of the simplified brake model.

	Component	Parameter	Value	Unit
		Thickness	25.	mm
		Inner radius	95.5	mm
	Disc	Outer radius Density	152.5 7200.	mm kg/m 3
		Young Modulus	130.	GP a
		Poisson's ratio (ν)	0.3
		Thickness	7.	mm
		Inner radius	96.5	mm
		Outer radius	151.5	mm
	Back plate	Angle	2π/7	radians
		Density	7300.	kg/m 3
		Young Modulus	170.	GP a
		Poisson's ratio (ν)	0.3
		Thickness	10.	mm
		Inner radius	96.5	mm
		Outer radius	151.5	mm
		Angle	2π/7	radians
		Slot width	5.	mm
	Lining	Density	2600.	kg/m 3
		Young Modulus (E x )	6.5	GP a
		Young Modulus (E z )	3.5	GP a
		Poisson's ratio (ν xy )	0.15
		Poisson's ratio (ν xz )	0.20
		Shear Modulus (G yz )	1.	GP a

Table 2

 2 Friction and Rayleigh coefficients for the brake system.

	Parameter

Table 4

 4 MAC and relative eigenvector errors for a selection of modes. Comparison results are given for three following homotopy order (H0-H1) to (H2-H3) of the three projection basis configurations.Modes No MAC (%) e ψ (%) MAC (%) e ψ (%) MAC (%) e ψ (%)

	i-k	nm=74 H0-H1	nm=74 H1-H2	nm=74 H2-H3
	46-46	91.	15.34	98.	6.90	100.	1.34
	47-47	85.	20.20	97.	8.81	99.	6.20
	48-48	0.	70.31	97.	8.81	99.	6.20
	49-49	4.	63.59	97.	9.36	100.	3.22
	50-50	91.	15.23	97.	9.36	100.	3.22
	51-51	99.	5.33	98.	8.57	97.	17.05
	52-52	95.	16.26	97.	14.84	94.	27.45
	53-53	94.	12.40	97.	9.31	100.	1.80
		nm=93 H0-H1	nm=93 H1-H2	nm=93 H2-H3
	46-46	97.	11.04	100.	4.02	100.	0.69
	47-47	1.	11.32	100.	9.30	100.	2.61
	48-48	72.	8.41	100.	8.14	100.	3.03
	49-49	0.	8.41	100.	8.14	100.	3.03
	50-50	97.	27.63	100.	42.86	100.	4.82
	51-51	97.	27.63	100.	42.86	100.	4.82
	52-52	95.	14.20	100.	5.89	100.	6.16
	53-53	92.	16.74	100.	5.89	100.	6.16
		nm=113 H0-H1	nm=113 H1-H2	nm=113 H2-H3
	46-46	97.	8.31	100.	2.27	100.	0.51
	47-47	1.	67.79	100.	1.03	100.	0.38
	48-48	72.	27.41	100.	1.48	100.	0.61
	49-49	0.	69.88	100.	1.48	100.	0.61
	50-50	97.	14.85	100.	3.15	100.	0.72
	51-51	97.	24.93	100.	3.15	100.	0.72
	52-52	95.	12.54	100.	1.51	100.	0.39
	53-53	92.	14.14	100.	1.51	100.	0.39

Acknowledgements

This work was carried out within the framework of the CNRS Research Federation on Ground Transports and Mobility, in association with the ELSAT2020 project supported by the European Community, the French Ministry of Higher Education and Research and the Hauts de France Regional Council. The authors gratefully acknowledge the support of these institutions.

Table 3

Real parts and natural frequencies for a selection of modes. Values are given for the three configurations of the projection basis from homotopy order 0 (H0) up to order 3 (H3).

Modes No

Re(s i )


François-Xavier Simon
email: francois-xavier.simon@inrap.fr

Mathias Pareilh-Peyrou

Solène Buvat

Alfredo Mayoral

Philippe Labazuy

Karim Kelfoun

Alain Tabbagh

Quantifying multiple electromagnetic properties in EMI surveys: a case study of hydromorphic soils in a volcanic contextthe Lac du Puy (France)

Keywords: EMI, VLF frequency range, conductivity, magnetic susceptibility and viscosity mapping, hydromorphic soils, igneous environment

Quantifying multiple electromagnetic properties in EMI

surveys: A case study of hydromorphic soils in a volcanic context -The Lac du Puy (France)

François-Xavier Simon, Mathias Pareilh-Peyrou, Solène Buvat, Alfredo

Mayoral, Philippe Labazuy, Karim Kelfoun, Alain Tabbagh

To cite this version:

François-Xavier Simon, Mathias Pareilh-Peyrou, Solène Buvat, Alfredo Mayoral, Philippe Labazuy, et al.. Quantifying multiple electromagnetic properties in EMI surveys: A case study of hydromorphic soils in a volcanic context -The Lac du Puy (France). Geoderma, 2020Geoderma, , 361, pp.114084. 10.1016Geoderma, /j.geoderma.2019.114084.114084. hal-02467986

Introduction

Palaeosols and sedimentary layers can provide valuable records of changing Holocene environmental conditions through careful multi-proxy geo-archaeological and paleoenvironmental analysis (e.g. [START_REF] Dotterweich | The history of soil erosion and fluvial deposits in small catchments of central Europe: Deciphering the long-term interaction between humans and the environment -A review[END_REF][START_REF] Dreibrodt | Are mid-latitude slopes sensitive to climatic oscillations? Implications from an Early Holocene sequence of slope deposits and buried soils from eastern Germany[END_REF][START_REF] Henkner | Archaeopedology and chronostratigraphy of colluvial deposits as a proxy for regional land use history (Baar, southwest Germany)[END_REF].

In this context, frequency domain electromagnetic induction measurements have proven to be a useful tool in the first stages of investigation and prior to coring or trenching [START_REF] Sudduth | Relating apparent electrical conductivity to soil properties across the north-central USA[END_REF][START_REF] Saey | Reconstructing the paleotopography beneath the loess cover with the aid of an electromagnetic induction sensor[END_REF][START_REF] Boaga | The use of FDEM in hydrogeophysics: A review[END_REF]. These techniques can provide information about electromagnetic properties and make it easier to identify the nature, geometry and thickness of superficial formations and sedimentary bodies containing paleo-environmental information.

Moreover, once drilling or trenching has been performed, geophysical observations calibrated with field data can be used to extrapolate stratigraphic or sedimentological observations [START_REF] Bendjoudi | Riparian wetlands of the middle reach of the Seine river (France): historical development, investigation and present hydrologic functioning, a case study[END_REF][START_REF] Hendrickx | Inversion of Soil Conductivity Profiles from Electromagnetic Induction Measurements: Theory and Experimental Verification[END_REF], Triantafilis and Monteiro Santos 2010[START_REF] Neely | Apparent electrical conductivity response to spatially variable vertisol properties[END_REF]).

In the framework of a broader geo-archaeological survey of a small pond, the Lac du Puy, situated at the top of the Puy de Corent in central France (Mayoral et al. 2018a), we used two EMI frequency domain instruments commonly used for environmental and soil studies.

First the EMP400 (GSSI Ltd, Nashua, NH, USA) where the inter-coil separation, L, equals 1.2 m, and secondly the EM31 (Geonics Ltd, Mississauga, Ontario, Canada) where the inter-coil separation equals 3.66 m, both operating with frequencies close to 10 kHz. Their geometry and operating frequency define the domain of electromagnetic induction used. Here it is the Very Low Frequency (VLF, 3-30 kHz) range, in which, regardless of the type of electromagnetic transmitted signal, the electrical conductivity, σ, is the dominant cause of the underground responses. Consequently, for both instruments, only the apparent electrical conductivity is usually considered [START_REF] Frischknecht | Profiling Methods Using Small Sources[END_REF]. However, four properties are actually involved (Fig. 1): two electrical (electrical conductivity and dielectric permittivity) and two magnetic (magnetic viscosity and magnetic susceptibility):

(1) While electrical conductivity represents the global motion of electric charges at the macroscopic scale, another property, permittivity, is required to describe the electrical polarization, which is the oscillations in electrical charge at a microscopic scale (from pore or clay platelet to molecular scales). In the absence of significant changes in salinity or temperature, variations in both properties reflect variations in the granularity, mainly the clay content, or in the soil water content. Both properties can be found in the Maxwell-Ampere equation in the expression,

r i    0 
, where ω is the angular frequency, ε0 the vacuum permittivity, εr the relative permittivity and 1   i . The relative weighting of the permittivity depends on the frequency: in the low frequency (EMI) domain it can only be strong if εr is strong [START_REF] Benech | Demonstrating the contribution of dielectric permittivity to the in-phase EMI response of soils: example of an archaeological site in Bahrain[END_REF][START_REF] Tabbagh | Permittivity mapping in the VLF-LF range using a multi-frequency EMI device: first tests in archaeological prospection[END_REF].

(2) As the acquisition of induced magnetization by a mineral grain is not immediate, there is a delay in acquisition or loss of the induced magnetization [START_REF] Dabas | Time domain magnetization of soils (VRM), experimental relationship to quadrature susceptibility[END_REF][START_REF] Mcneill | The magnetic susceptibility of soils is definitely complex[END_REF][START_REF] Jordanova | Soil Magnetism. Applications in Pedology[END_REF]). This delay is taken into account by considering that the magnetic susceptibility is a complex quantity,

qu ph i    
, where κph is the inphase magnetic susceptibility (hereafter simply called magnetic susceptibility) and κqu is the magnetic viscosity. The magnetic viscosity depends on the size of the magnetic grains. It is strong for small single domain magnetite grains (in the [10 -30 nm] equivalent diameter range), small for larger single domain grains (around 1 µm) and once again strong for large multi-domain grains (20 µm and greater). These magnetic properties can be used jointly to describe the degree of pedogenesis, as the quantity and size of magnetic grains are directly linked to redox processes that affect iron oxides/hydroxides. Although the presence of organic matter affects these processes [START_REF] Garcia | Soil geochemical methods in archaeo-geophysics: Exploring a combined approach at sites in Scotland[END_REF], the content of iron oxides of the parent rock plays a key role in the magnitude of soil magnetic susceptibility. In general, high magnetic susceptibility would be expected for basaltic rocks (1000 to10000 10 -5 SI).

Frequency domain EMI apparatuses with a coil spacing of greater than 3 m and a frequency lower than 30 kHz -such as the EM31-have been used in a wide variety of environmental studies since 1979 [START_REF] De | Measurement of Apparent Electrical Conductivity of Soils by an Electromagnetic Induction Probe to Aid Salinity Surveys[END_REF]. Furthermore, over 200 papers have been published using the conductivity measurements obtained with this device [START_REF] Doolittle | The use of electromagnetic induction techniques in soils studies[END_REF], but to the best of our knowledge only very few studies look at the other EM property measurements for this instrument. On the other hand, in the field of archaeological prospection, where the separation between the coils remains limited, less than two meters, it has been noted that the magnetic properties play an important part in ground response (Tite and Mullins 1970[START_REF] Scollar | Archaeological prospecting and remote sensing[END_REF][START_REF] Sarris | Mapping of quadrature magnetic susceptibility/magnetic viscosity of soils by using multi-frequency EMI[END_REF][START_REF] Smedt | Comparing apparent magnetic susceptiblity measurements of a multi-receiver EMI sensor with topsoil and profile magnetic susceptibility data over weak magnetic anomalies[END_REF].

In the Puy de Corent pond, the presence of basaltic rocks (basanite) underlying the depression is of particular interest in the study of magnetic properties and electrical polarization phenomena. Such rock types usually contain a high content of stable single domain magnetic grains with high susceptibility [START_REF] Grison | Magnetic and geochemical characterization of Andosols developed on basalts in the Massif Central, France[END_REF], very low viscosity and significant electrical polarizability. Therefore, the purpose of our work was to investigate the four measurable electromagnetic properties rather than using electrical conductivity alone, and to assess their relative abilities to improve interpretation of three-dimensional sedimentary infill (Triantafilis et al., 2013a, Triantafilis et al., 2013b), well documented in Mayoral et al. 2018a but mainly based on direct and spatially limited stratigraphic observations. This non-intrusive approach provides benefits in terms of both time and resources. The raw data is collected by the two instruments and then inverted using point by point 1D modelling, considering first the electrical conductivity alone, and then incorporating the electrical conductivity, the complex magnetic susceptibility and the relative dielectric permittivity in a second round of modelling.

Material and methods

2-1. Study area and Stratigraphy

The study site is located on top of the Puy de Corent (621 m.a.s.l.), a volcanic plateau dominating the calcareous Limagne Plain in central France. The plateau has been shaped by erosion and topographic inversion of a 3 million year-old scoria cone and an associated basaltic lava flow formed during the Pliocene volcanism of the Massif Central [START_REF] Greffier | Aspects géomorphologiques et stabilité des versants au sud de Clermont-Ferrand[END_REF][START_REF] Nehlig | Les volcans du Massif central[END_REF]. In the quasi-flat lower part of the plateau, a small circular depression of circa. 2 ha forms a wetland representing a drained pond (Fig. 2). This depression was formerly interpreted as being a secondary crater, but is more likely a pseudo-sinkhole, caused by subsidence phenomena affecting the basaltic lava flow, likely triggered by dissolution of underlying Oligocene sedimentary rocks including marls, limestones and gypsum (Mayoral et al., 2018). Litho-stratigraphic features of the Lac du Puy have already been described in detail in previous work (Mayoral et al. 2018), and are summarized here in Figure 2. The bedrock is basaltic and is formed of the Puy de Corent lava flow. The sedimentary in-fill of the basin is mainly clayey and can reach a thickness of up to 1.8 or 2 m. The lower stratigraphic units have a clay matrix which is very rich in coarse volcanic material (sand to gravel), resulting from detrital input from the plateau and local bedrock weathering and disaggregation. The middle stratigraphic units are dominated by massive, heavy clays which have undergone various phases of pedogenesis since the Neolithic, with dominant hydromorphic features (gleyic and vertic features) and rare coarse particles. The upper units are mainly clayey loams including several basaltic fragments and pottery sherds, and show Roman to present-day soil development.

Magnetic susceptibility was measured on two different cores along a trench and is summarized as an average per sedimentary unit. Susceptibility values were measured with the single frequency MS2-E (Bartington Ltd, UK, 1cm spatial resolution, 2 kHz) on the fine soil fraction. This resulted in relatively low susceptibility values (less than 50.10-5 SI), not taking into account the effect of the basaltic and sherd fragments.

2-2. Devices and survey conditions

Two different EMI devices were used on the study site. The EMP400 operates at three frequencies below 16 kHz (5, 8 and 15 kHz were selected). This instrument has one transmitting coil and one receiving coil at a distance of L=1.21 m from the transmitting coil. The coils are coplanar, and measurements can be made in either the HCP (horizontal coplanar) or the VCP (vertical coplanar), which is configured by rotating the instrument itself. The depths of investigation for the EMP400 differ between the HCP and VCP, and for the electrical and magnetic properties (Tabbagh 1986b). For the electrical properties, it is generally 0.8L in VCP and 1.5L in HCP, whereas for the magnetic ones the distances are 0.6L in VCP and 0.2L in HCP, with a change in the sign of the response for a clearance h=0.38 L.

The second device, the EM31, has a coplanar configuration, a 3.66 m inter-coil separation and a single frequency of 9.8 kHz. The operator can hold the instrument at shoulder height with h=1 m clearance above the ground surface. It can also be pulled using a trolley at a lower clearance. The depths of investigation follow the same general guidelines for L as for the EMP400.

The Lac du Puy area was surveyed with the two instruments in both HCP and VCP configurations, in combination with a GNSS positioning system with a sub-decimetric accuracy in order to locate the instrument operator precisely. A profile was collected every two meters with a frequency of data collection of 1 Hz for the EM31 and 10 Hz for the EMP400. Data for both instruments and configurations were collected on October 2015, resulting in eight different independent sets of data: EM31 in HCP, EM31 in VCP, EMP400 in HCP for 5, 8 and 15 kHz and EMP400 in VCP for 5, 8 and 15 kHz. Each was resampled with the same 2 m x 2 m grid meshing in order to enable simultaneous interpretations and inversions using these different sets.

2-3. 1D forward modelling

We present here the 1D forward modelling, which takes the four electromagnetic properties into account (see also Ward and Homann 1988). It allows us to define the apparent properties when the ground is assumed to be homogeneous, and to inverse for properties and thicknesses of the successive layers when a multi-layer model is considered.

 For horizontal coplanar (HCP), where the transmitter is a vertical magnetic dipole, Mz, located at (0, 0, -h) and the receiver a small coil located at (0, L, -h),

    0 0 2 2 ) ( ) ( 4 0      d L J e R M H h u v z z (1)
 For vertical coplanar (VCP), where the transmitter is a horizontal magnetic dipole, Mx, located at (0, 0, -h) and the receiver a small coil located at (0, L, -h),

               0 0 1 2 1 0 2 2 0 ) ( ) ( 1 ) ( ) ( 4 0 0          d L J e R L d L J e R M H h u h u H x x (2)
In these expressions J0 and J1 are the Bessel's function of the first kind,

2 0 2 0     u with 2 0 0 2 0      
, and

1 0 0 1 0 0 0 ) ( Y u Y u u R v       
, Y1 being recursively calculated by starting at the top of the deepest layer, N, by

N N N u Y   
and using the formula:

i i i i i i i i i i i i i e u Y u e u u Y u Y tanh tanh 1 1        
. In this formula, ei is the thickness of the i th layer, and

2 2 i i u     with 2 2        i i i i i i   and 1   i . ) ( ) ( ) ( 2 0 2 0 1     v H R u R R   . The function 𝑅 𝐻 (𝜆) is calculated by 1 0 0 1 0 0 0 ) ( Z i u Z i u u R H        , Z1
being recursively calculated by starting at the top of the deepest layer by

N N N N i u Z     
and using the formula:

i i i i i i i i i i i i i i i i e u Z i u e u i u Z i u Z tanh tanh 1 1              
. To get results which are independent of the transmitted power, the secondary field expressions are usually normalized by the primary amplitude, Hp=M/(4πL 3

) and the measurements are expressed in terms of Hs/Hp , in ppm (parts per million) or in ppt (parts per thousand) depending on the device manufacturer (ppm for the EMP400 and ppt for the EM31, 1ppt=1000 ppm).

To illustrate the physical meaning of these complex analytic expressions by showing the influence of any property or instrument parameter, approximations (Thiesson et al. 2014) can be calculated based on the range of values of the different properties. The static approximation ( 0   ) establishes that the ratio of the secondary field to the primary field for a coplanar instrument (HCP or VCP) located at h=0 clearance above the surface would be

2 qu ph p s i H H     
. This suggests that the magnetic susceptibility response would be in-phase, the viscosity response in quadrature and their dependence on frequency and coil separation limited. Furthermore, the Low Induction Number (LIN) approximation,

1 2  L  , which corresponds, if µ=µ0 and h=0, to 4 2 0 L i H H p s    
suggests that the conductivity effects would mainly be in quadrature and proportional to the conductivity, the frequency and the square of the coil separation (this approximation is the one commonly used by most instrument manufacturers to define ECa).

If we take into account the permittivity, this last approximation becomes 4 4

2 2 0 0 2 0 L L i H H r p s         
, which suggests that the permittivity response would be first order in-phase and would drastically increase with the frequency and coil separation.

Consequently, the in-phase response can be dependent on both the magnetic susceptibility and the relative permittivity. Note also that even if the displacement currents are negligible in comparison with conduction ones,

    r i 0 
, the permittivity response may be significant in comparison with the susceptibility response.

2-4. Instrument calibration and sensitivity analysis

At each geometrical configuration and frequency an instrument measures two quantities, the in-phase and quadrature responses, and it is important to understand the influence of each of the four properties (electrical conductivity (or Eca), dielectric permittivity, magnetic susceptibility (or Msa) and magnetic viscosity) on these responses.

To assess this, the calculated responses, as a function of each of the four properties, are expressed by the ratios of the secondary field to the secondary field of a reference ground, where σ=0.01 Sm -1 , εr=200, κph=100 10 -5 SI and κqu =0.06κph. In this trial each property is varied from 0.1 to 10 times that of the value used in the reference case value. Figure 3 presents the in-phase modelled responses for the EMP400, and Figure 4 those for the quadrature situation for both configurations and each frequency. For the low-frequency instrument, the influence of permittivity is negligible, as illustrated by the flat curve in Figure 3 at all configurations and frequencies. The magnetic susceptibility dominates the in-phase response with a linear dependence (Fig. 3), and conductivity controls the quadrature responses with a quasi-linear dependence (Fig. 4). These characteristics are in complete agreement with all that has already been published [START_REF] Mcneill | Electromagnetic terrain conductivity measurement at low induction numbers[END_REF], Tabbagh 1986a). However, a precise interpretation needs to take the effect of the conductivity on the in-phase response into consideration, which increases, but not in direct proportion to the conductivity (Fig. 3), and the effect of the magnetic viscosity on the quadrature response (Fig. 4). The procedure to determine these properties is presented in [START_REF] Sarris | Mapping of quadrature magnetic susceptibility/magnetic viscosity of soils by using multi-frequency EMI[END_REF].

No drift has been observed for the EMP400 (but the temperature changes during the acquisition remained limited). A major drawback of the raw data sets is that they first require the determination of the zero offsets or signal shifts. For in-phase components such shifts are mainly caused by mechanical deformations that modify the direct influence of the primary field, but, surprisingly, offsets also exist in quadrature. For each configuration, this step was achieved by a trial and error process to overcome this problem. In quadrature, the different apparent conductivity values, calculated at each frequency and by considering the differences between each frequency pair, must be as close as possible to each other and positive, while the magnetic viscosity must be positive [START_REF] Sarris | Mapping of quadrature magnetic susceptibility/magnetic viscosity of soils by using multi-frequency EMI[END_REF]. In-phase, the different susceptibility values must be as close as possible and positive. The zero offsets are given in Table 1.

Following the same approach as for the EMP400, Figures 5 and6 show the evolution of the ratio of the in-phase and quadrature secondary fields to those of the reference case for the EM31, in which each property varies from 0.1 to 10 times the value used in the reference case. The dependence of the quadrature response on magnetic viscosity is, in this case (h=1m), quasi null (Fig. 6); therefore, the conductivity determination can be directly achieved with this component at a single frequency. The in-phase component is sensitive to conductivity, permittivity and inphase magnetic susceptibility (Fig. 5), the permittivity response is illustrated by a quasi-flat curve for all configurations and clearance; it is significantly smaller than those of the susceptibility and conductivity. Determination of the apparent properties thus follows a twostep approach: (1) determination of the apparent conductivity using the quadrature response;

(2) taking into account the in-phase electrical conductivity response, determination either of the apparent magnetic susceptibility or of both the permittivity and the susceptibility if both HCP and VCP configurations are available. This last step follows the linear process described in [START_REF] Benech | Demonstrating the contribution of dielectric permittivity to the in-phase EMI response of soils: example of an archaeological site in Bahrain[END_REF].

With EM31, no zero offset correction is necessary for the quadrature part of the signal because of the good field calibration of the device achieved by the manufacturer and of the large coil spacing which increases the signal/noise ratio. On the other hand, the zero offset corrections for the in-phase part of the signal must be evaluated (this component is used in metal detection). Here, following the same the calibration protocol as for the EMP400, in-phase signal offsets were established at 56 ppt for HCP (56000 ppm) and 64 ppt (64000 ppm) for VCP by a trial and error process. These high values correspond to a significant uncertainty regarding the exact offset value, but not for the spatial variations of the properties.

3-Results

3-1. Apparent property maps

Each apparent property is defined here as the property of homogeneous flat ground that would give the same response as that obtained with the considered instrument.

The maps resulting from the EMP400 measurements are presented in Figure 7 for the apparent conductivity, and in Figure 8 for the apparent magnetic susceptibility. As the magnetic viscosity, the conductivity and the depth of investigation are independent of the frequency, the difference in quadrature responses between the 15 kHz and 5 kHz frequencies allowed us to remove the magnetic viscosity effect [START_REF] Sarris | Mapping of quadrature magnetic susceptibility/magnetic viscosity of soils by using multi-frequency EMI[END_REF]. To determine the electrical conductivity we thus compared these differences to the corresponding theoretical reference curves calculated from equations (1) for HCP, or (2) for VCP. We also used theoretical reference curves for each frequency to determine the electrical conductivity using the quadrature responses. Using this electrical conductivity value, it is then possible to extract the magnetic viscosity value, which is shown for the VCP and HCP configurations in Figure 9.

As illustrated in Figure 7, the maps of electrical conductivity at the different frequencies show identical results, in complete agreement with the values of the Induction Number that establish the similarity in depths of investigation. This result corroborates the solution adopted for the computation of magnetic viscosity. The HCP and VCP maps for apparent conductivity are also very similar. The HCP configuration exhibits a greater range of values, and lower values at the edge of the depression, in agreement with an investigation at greater depth. In the center, the values are higher than 0.04 S.m -1 and, as expected, confirm the clayey nature of the filling and clearly delineate the topographic trough.

For magnetic properties, reference curves deduced from equations ( 1) and ( 2) are used to determine the apparent magnetic susceptibility, taking into account the estimated electrical conductivity at each point. The HCP and VCP susceptibilities are globally very high, between 300 and 3,800 10 -5 SI, in agreement with the presence of a basaltic rock substratum (Trigui and Tabbagh 1990). The differences between results for the three frequencies are small, as predicted by the theory. The strong similarity confirms the efficiency of the correction of the induction effect on the in-phase signal for HCP and VCP configurations, respectively. However, the sediment in-fill has lower susceptibility values than the border, at around 700 10 -5 SI. The horizontal delineation is in complete agreement with that of the conductivity. The magnetic viscosity is relatively low and its lateral distribution differs from that of the susceptibility. In HCP the magnetic viscosity (depth of investigation of 0.2 L) is higher than in VCP (depth of investigation: 0.6 L) which suggests that it is stronger for the superficial layer of the depression. This observation is in complete agreement with the field observations. The uppermost 45-50 centimeters, of slightly varying thickness, contain many sherds and are very different to the deeper layer (Mayoral et al. 2018a). The ratio ph qu   values, being higher than 4% inside the topographic trough and lower than 1% at the border, confirm that the basalt mainly contains stable single-domain magnetic grains (high susceptibility and low magnetic viscosity) while the in-fill has been enriched by pedogenic processes -such as redoximorphic processes or alteration of mineral grains (Mayoral et al. 2018) -into smaller single-domain grains closer in volume to the superparamagnetic range [START_REF] Maher | Formation of ultrafine-grained magnetite in soils[END_REF]. By comparison, a ratio of 6% or more is expected for common non-hydromorphic soils [START_REF] Dabas | Time domain magnetization of soils (VRM), experimental relationship to quadrature susceptibility[END_REF]. As only a single frequency device (MS2-E) was used on the laboratory soil sample measurements, no direct comparison in terms of frequency variation for the susceptibility can be carried out for this field observation.

HCP and VCP conductivity maps obtained from the EM31 quadrature signal (Figure 10) confirm the horizontal distribution shown by the EMP400 maps. As expected, given that the investigation depth is greater, the conductivity values are about one and a half times lower than for the EMP400.

For the in-phase signal, due to the longer coil separation, the effects of both the magnetic susceptibility and the relative permittivity should be considered. The raw data variability is presented in Table 2. The sensitivities of EM31 to the homogeneous ground properties are given in Table 3. The lower sensitivity to κph of the HCP response is due to the instrument's clearance, h=1 m, because this elevation is close to the change in sign which occurs in the HCP configuration for magnetic responses (at h=1.39 m for a coil separation of 3.66 m). If we calculate the apparent magnetic susceptibility from the in-phase responses without considering the relative permittivity, the discrepancy between HCP and VCP results is marked (Figure 11):

the spatial distributions are significantly different with a higher dynamic in HCP (in the range 7,500 to 50,000 10 -5 SI), and a lower dynamic in VCP (in the range 2,500 to 18,000 10 -5 SI).

As potentially imperfect offset corrections would increase this effect, the significant difference between the two configurations can probably be explained by the sensitivity to susceptibility distribution in relation to the depth.

The maps presented in Figure 12 show the results with both the permittivity and the susceptibility taken into account to define the homogeneous half-space properties using both HCP and VCP in-phase values [START_REF] Benech | Demonstrating the contribution of dielectric permittivity to the in-phase EMI response of soils: example of an archaeological site in Bahrain[END_REF]. In these apparent property maps, the susceptibility variations remain similar to those observed in Figure 11, whereas the permittivity variations are huge, with a 140,000 interquartile deviation and a big gap for the zero offset. The range of values looks unlikely and it can only be assumed that the basalt border has lower values than the sedimentary in-fill in the depression. The high susceptibility values agree with the presence of a basaltic substratum underlying the sedimentary in-fill and are logically greater than the values delivered by the EMP400. The permittivity values logically result from the low sensitivity presented in Table 3. Considering the high magnetic context and the unlikely apparent permittivity values, the intervention of a permittivity response looks unrealistic. This conclusion will be reconsidered later on when considering a layered terrain.

From the apparent property mapping, it can be seen that the filling is more conductive and possibly more polarizable than the surrounding basalt. Its conductivity indicates that it is dominated by fine particles. It has a lower magnetic susceptibility with respect to the surrounding weathered basalt, but a higher magnetic viscosity. Taking all this into account, the electromagnetic properties can provide an accurate description of magnetic behavior in addition to the electrical properties. We are then able to describe not only the clay filling of the depression but also some pedogenetic phenomena, and even provide further details about the lateral and horizontal variations in the hydromorphic soils in a volcanic context, which is poorly revealed by a single electrical conductivity mapping based approach.

3-2 1D inversion processing using quadrature data

Based on the 1D forward modelling [START_REF] Guérin | Interpretation of Slingram conductivity mapping in near-surface geophysics: using a single parameter fitting with 1D model[END_REF], we carried out 1D inversions, firstly considering the electrical conductivity alone, then taking all the electromagnetic properties of the soil into account. Due to the quality of the data, it was not necessary to select a lateral smoothing inversion process, as it was possible to achieve this point by point. The forward operator, G, uses equation ( 1) for HCP data, and equation ( 2) for VCP data. At each point, an iterative damped least square linearization was applied, which started from the m0 vector of the a priori values of the conductivity and/or thickness and calculated Δm, the solution of the equation,

) ( ) ( 1      i T T Gm d J m I J J 
, mi-1 being the vector solution at the (i-1) th iteration, d the vector of the data, I the identity matrix, J the Jacobian of G, and λ the regularization parameter. λ has a starting value equal to half of the trace of the J T J matrix.

To carry out a more in depth characterization of the sedimentary in-fill, we attempted to reconstruct the depth to the bedrock using a 1D inversion in a two-step approach. In the first step, we only considered the HCP and VCP quadrature responses by using the differences between 15 kHz and 5 kHz for the EMP400 quadrature data, and the quadrature data at 9.8 kHz for the EM31. In the second step, we considered the four properties using eight independent sets of data: quadrature and in-phase. For the conductivity profile, following comparison with the field data (Mayoral et al. 2018), we used a three-layer model (see Figure 2) with a superficial soil layer (the two top clay loam layers), a second layer corresponding to the sedimentary infill (massive heavy clay layer), and a third layer corresponding to the possibly weathered basaltic rocks (sandy bottom and weathered basaltic rock layers). As an initial model of the parameters, the superficial layer had a conductivity of σ1=20 mS.m -1 (or 50 Ω.m resistivity) and a thickness of e1=0.15 m; the sedimentary in-fill a conductivity of σ2=83.3 mS.m -1 (or 12 Ω.m resistivity) (a little bit higher than the EMP400 minimum) and a thickness of e2=1 m; and the basement a conductivity of σ3=3.33 mS.m -1 (or 300 Ω.m resistivity) (lower than any of the EM31 HCP measurements). The conductivities and thicknesses of the different layers were tested and it was observed that the parameters of the uppermost layer had only a very limited influence. EMI inversions were also subject to so-called 'cases of equivalence': in the case of a conducting layer between two more resistive ones (H type for vertical electrical soundings), the conductive layer can be exchanged for any other layer with the same  3 2 e ratio [START_REF] Guérin | Interpretation of Slingram conductivity mapping in near-surface geophysics: using a single parameter fitting with 1D model[END_REF]. Finally, by allowing only e2 and σ1 to vary, and keeping the other parameters fixed, the thickness distribution presented in Figure 13 was obtained.

These results suggest that the 0.75 m thickness value corresponds to the contour of the depression itself, due to a thickness contrast with the surrounding second layer. This thickness variation is in complete agreement with the profile published by [START_REF] Guillemoteau | 1D sequential inversion of portable multi-configuration electromagnetic induction data[END_REF] when interpreting DualEM421 results and with the field observation (Mayoral et al. 2018). We also carried out a control of the inversion results by comparing field observations with the results of the inversion using the quadrature part of the signal. Figure 14(a) shows that the estimations of the thicknesses of the three layers made it possible to clearly describe the depth of the high clay content layers filling the depression. Only one of the logs (log 3) located in the central trench of the site (Figure 14), does not fit the inversion results. This difference could derive from the 1D interpretation which leads to a smoothed model of the bottom of the second layer where there is in reality a local rise in the basalt.

3-3 1D inversion using both quadrature and in-phase data sets

Using the three-layer model obtained by 1D modelling of the conductivity results, the other EM properties for each layer were determined point by point. The same issue as above recurs: is it possible to determine both the susceptibility and permittivity or can the influence of the permittivity be neglected? In other words, which hypothesis gives the more realistic results with respect to the geological/pedological context and the external information available describing the depression as having an irregular base shaped into several alteration zones?

When the permittivity is neglected, the magnetic susceptibility obtained for the second layer by using four sets of in-phase data is presented in Figure 15. The permittivity for the three layers is fixed at an ordinary value of: εr1=εr2=εr3=200, and with κph1=500 10 -5 SI and κph3=2000

10 -5 SI. The map is in agreement with the apparent susceptibility maps in terms of the lateral distribution and the parameter's magnitude (Figures 8 and 11).

Figure 16 presents the results for the case where the permittivity of the second layer is variable. These results are obtained with εr1=εr3=2000 but are practically independent of the first and third layer fixed permittivity values. The susceptibility is reduced by a factor close to two and, again, the permittivity values are very high and unrealistic. For a permittivity of 10 6

(median value obtained), the quantity ωεrε0 would equal 0.54 Sm -1 , seven times greater than the conductivity itself (83 mS.m -1 ). Considering this extremely high and unlikely estimated range of values of dielectric permittivity, the second hypothesis, where the variations of the permittivity are not characterized for the EMP400 and EM31 instruments, and where the inphase results correspond to magnetic susceptibility only, must be adopted.

A comparison with the susceptibility observed at a scale of 1 cm in the cores ('U-Channel' cores), which is summarized in Figure 14, shows a significant discrepancy between the magnitude of inverted in-field susceptibility values and the laboratory ones (MS2-E) for the in-fill of the depression. This discrepancy is explained if one considers that the presence of basalt gravels and pebbles (with around 12000 10 -5 SI susceptibility) is taken into account by in-field measurements rather than laboratory one on the cores. These coarse particles would correspond to one fourth or one fifth of the sediment volume at a metric scale.

Conclusion

This case study confirms the relevance of this type of light EMI apparatus for the rapid exploration of superficial sedimentary/pedological formations, and will increase the reliability of future sediment sampling at the Lac du Puy. Two different single receiver instruments were used, although a multi-receiver/multi-frequency apparatus would produce the same results and be easier to use.

The highly magnetic context due to the basaltic substrate led us to look at the in-phase measurements and try to describe the magnetic susceptibility and magnetic viscosity of the infill. This process is more common for shorter separation instruments, but nothing has yet been published about this usage for the EM31. The superficial soil and the depression in-fill have a lower magnetic susceptibility and a higher magnetic viscosity than the surrounding weathered basalt and their spatial distribution is slightly different to that of the sediment thickness. Due to the improbable magnitude of the inverted permittivity results, no role can be attributed to the dielectric permittivity variations and the full signal can be attributed to the complex magnetic susceptibility and electrical conductivity. Further studies should be done in the future to assess the effect of permittivity on the in-phase EM31 signal. In this volcanic substrate case study, magnetic susceptibility responses outweigh the electrical polarization when using the EM31. It would be worth doing further investigations in order to assess the potential for this type of instrument in other environments, such as coastal areas or salt-rich soils.

In conjunction with the field characterization of the sedimentary accumulation, measurements of the quadrature and in-phase components of the Hs/Hp ratio using the HCP and VCP configurations for both instruments provided a clear description of the conductivity variations. It enabled the Lac du Puy sedimentary in-fill to be identified, by mapping the shape and morphology of the bedrock, lending support to the hypothesis that the depression is a pseudo-sinkhole (Mayoral et al 2018a).

The analysis of a combination of magnetic properties allowed the 3D morphology of the sedimentary in-fill to be determined across the basin, and provided stratigraphic and sedimentological information. In addition, selected magnetic properties provided information on the spatial variation of some key features related to hydromorphic processes, in particular clayey granularity and the development of iron oxides/hydroxides. These findings, which will be further tested and refined in other parental rocks, could provide new means of rapid spatial and pedological characterization of hydromorphic soils and palaeosoils, with useful applications for palaeoenvironmental studies. Trigui M., [START_REF] Scollar | Archaeological prospecting and remote sensing[END_REF]. Magnetic susceptibilities of oceanic basalts in alternative fields.

Journal of Geomagnetism and Geoelectricity, 42, 621-636.

Ward S. H., Homann G. W., 1988. Electromagnetic Methods in Applied Geophysics: Volume 1, Theory, M. N. Nabighian ed., SEG Publications, Tulsa OK. Figure 16: Inversion of both the magnetic susceptibility and the dielectric permittivity of the second layer.

Figure captions

Tables :

Table 1: EMP400 zero offset corrections in ppm

Figure 1 :

 1 Figure 1: Schematic view of the relationship between the soil properties, the geophysical

Figure 2 :

 2 Figure 2: a) Site location in Central France, b) Aerial view of the Plateau de Corent and the Lac

Figure 3 :

 3 Figure 3: EMP400 (L=1.21 m, h=0.25 m) ratio of the in-phase secondary field to the in-phase

Figure 4 :

 4 Figure 4: EMP400 (L=1.21 m, h=0.25 m) ratio of the quadrature secondary field to the

Figure 5 :

 5 Figure 5: EM31 (L=3.66 m, f=9.8 kHz) ratio of the in-phase secondary field to the in-phase

Figure 6 :

 6 Figure 6: EM31 (L=3.66 m, f=9.8 kHz) ratio of the quadrature secondary field to the quadrature

Figure 7 :

 7 Figure 7: EMP400 apparent conductivity maps: determined using the difference in quadrature

Figure 8 :

 8 Figure 8: EMP400 apparent magnetic susceptibility maps, κph, at the three frequencies used (5,

Figure 9 :

 9 Figure 9: EMP400 apparent magnetic viscosity map κqu for HCP (left) and VCP (right)

Figure 10 :

 10 Figure 10: EM31 apparent electrical conductivity maps for HCP (left) and VCP (right)

Figure 11 :

 11 Figure 11: EM31 apparent magnetic susceptibility, κph, maps determined without taking the

Figure 12 :

 12 Figure 12: EM31 maps of the apparent magnetic susceptibility, κph, and of the apparent relative

Figure 13 :

 13 Figure 13: Thickness of the second layer, e2, resulting from a 1D inversion when σ1 and e2 are

Figure 14 :

 14 Figure 14: a) comparison of field observations and results based on the estimation of the

Figure 15 :

 15 Figure 15: Inversion of the magnetic susceptibility of the second layer, ignoring permittivity

Figure 1 :

 1 Figure 1: Schematic view of the relationship between the soil properties, the geophysical

Figure 2 :

 2 Figure 2: a) Site location in Central France, b) Aerial view of the Plateau de Corent and the Lac

Figure 3 :

 3 Figure 3: EMP400 (L=1.21 m, h=0.25 m) ratio of the in-phase secondary field to the in-phase

Figure 4 :

 4 Figure 4: EMP400 (L=1.21 m, h=0.25 m) ratio of the quadrature secondary field to the

Figure 5 :

 5 Figure 5: EM31 (L=3.66 m, f=9.8 kHz) ratio of the in-phase secondary field to the in-phase

Figure 6 :

 6 Figure 6: EM31 (L=3.66 m, f=9.8 kHz) ratio of the quadrature secondary field to the quadrature

Figure 7 :

 7 Figure 7: EMP400 apparent conductivity maps: determined using the difference in quadrature

Figure 8 :

 8 Figure 8: EMP400 apparent magnetic susceptibility maps, κph, at the three frequencies used (5,

Figure 9 :

 9 Figure 9: EMP400 apparent magnetic viscosity map κqu for HCP (left) and VCP (right)

Figure 10 :

 10 Figure 10: EM31 apparent electrical conductivity maps for HCP (left) and VCP (right)

Figure 11 :

 11 Figure 11: EM31 apparent magnetic susceptibility, κph, maps determined without taking the

Figure 12 :

 12 Figure 12: EM31 maps of the apparent magnetic susceptibility, κph, and of the apparent relative

Figure 13 :

 13 Figure 13: Thickness of the second layer, e2, resulting from a 1D inversion when σ1 and e2 are

Figure 14 :

 14 Figure 14: a) comparison of field observations and results based on the estimation of the

Figure 15 :

 15 Figure 15: Inversion of the magnetic susceptibility of the second layer, ignoring permittivity

Figure 16 :

 16 Figure 16: Inversion of both the magnetic susceptibility and the dielectric permittivity of the

  Thiesson J.,Kessouri P., Schamper C., Tabbagh A. 2014. Calibration of frequency-domain electromagnetic devices used in near-surface surveying. Near SurfaceGeophysics, 12,[481][482][483][484][485][486][487][488][489][490][491]

	Tite M. S., Mullins C. E. 1970. Electromagnetic prospecting on archaeological sites using a soil
	conductivity meter. Archaeometry, 12, 97-104.
	Triantafilis J., Monteiro Santos F. A., 2010. Resolving the spatial distribution of the true
	electrical conductivity with depth using EM38 and EM31 signal data and a laterally constrained
	inversion model. Australian Journal of Soil Research, 48, 434-446.
	Triantafilis, J., Terhune IV, C.H., Monteiro Santos F.A., 2013a. An inversion approach to
	generate electromagnetic conductivity images from signal data. Environmental Modelling and
	Software, 43,88-95.
	Triantafilis J., Monterio Santos, F.A., 2013b. Electromagnetic conductivity imaging (EMCI) of
	soil using a DUALEM-421 and inversion modelling software (EM4Soil). Geoderma, 211-212

Table 2 :

 2 Variabilities of the in-phase raw data from EM31

Table 3 :

 3 Sensitivities of the EM31 in-phase responses to homogeneous ground susceptibility and permittivity when h=1 m

Table 2 :

 2 Variabilities of the in-phase EM31 raw data

EM31 (h=1m) HCP VCP CHE= 0.110 ppm CVE= 0.077 ppm CHK= -1.044 ppm/1 10 -5 SI CVK= 3.36 ppm/1 10 -5 SI

Table 3 :

 3 Sensitivities of the EM31 in-phase responses to a homogeneous ground susceptibility and permittivity when h=1 m