Thiesson

3D linear inversion of magnetic susceptibility data acquired by frequency domain EMI

Keywords: Magnetic susceptibility of soils, frequency domain EMI, 3D inversion of inphase susceptibility measurements

Low induction number EMI instruments are able to simultaneously measure a soil's apparent magnetic susceptibility and electrical conductivity. This family of dual measurement instruments is highly useful for the analysis of soils and archaeological sites. However, the electromagnetic properties of soils are found to vary over considerably different ranges:

whereas their electrical conductivity varies from ≤ 0.1 to ≥ 100 mSm -1 , their relative magnetic permeability remains within a very small range, between 1.0001 and 1.01 SI. Consequently, although apparent conductivity measurements need to be inverted using non-linear processes, the variations of the apparent magnetic susceptibility can be approximated through the use of linear processes, as in the case of the magnetic prospection technique.

Our proposed 3D inversion algorithm starts from apparent susceptibility data sets, acquired using different instruments over a given area. A reference vertical profile is defined by considering the mode of the vertical distributions of both the electrical resistivity and of the magnetic susceptibility. At each point of the mapped area, the reference vertical profile response is subtracted to obtain the apparent susceptibility variation dataset. A 2D horizontal Fourier transform is applied to these variation datasets and to the dipole (impulse) response of each instrument, a (vertical) 1D inversion is performed at each point in the spectral domain, and finally the resulting dataset is inverse transformed to restore the apparent 3D susceptibility variations.

It has been shown that when applied to synthetic results, this method is able to correct the apparent deformations of a buried object resulting from the geometry of the instrument, and to restore reliable quantitative susceptibility contrasts. It also allows the thin layer solution, similar to that used in magnetic prospection, to be implemented. When applied to field data it initially delivers a level of contrast comparable to that obtained with a non-linear 3D inversion. Over four different sites, this method is able to produce, following an acceptably short computation time, realistic values for the lateral and vertical variations in susceptibility, which are significantly different to those given by a point-by-point 1D inversion.

Introduction

Soil is produced by various complex processes, to which human activities can make a significant contribution. A complete, continuous and non-invasive description of a soil's structure is thus of primary importance in terms of improving our knowledge of ancient societies, and ensuring more relevant management of the current environment. Geophysical surveys have thus been recognized as an indispensable tool for subsurface environmental studies [START_REF] Butler | Near Surface Geophysics (Investigations in Geophysics)[END_REF][START_REF] Viscarra Rossel | Proximal soil sensing[END_REF]). In addition to its electrical resistivity, the soil's magnetic properties have sufficient variability, and are sufficiently related to past and present active pedological processes (namely redox) [START_REF] Evans | Environmental Magnetism, Principles and Applications of Enviromagnetics[END_REF]Heller 2003, Liu et al. 2012), to motivate their utilization in ground prospection campaigns. Magnetic prospection has been the most commonly used technique in archeological applications [START_REF] Aspinall | Magnetometry for archaeologists[END_REF], in which spatial variations of the Earth's magnetic field, and thus all variations in the ground's total magnetization, are measured. However, this technique fails to describe the vertical layering of a terrain, and is thus less useful for pedological applications in which it is important to identify the magnetic properties of each horizon; moreover, its ability to detect lens-like features is very poor [START_REF] Scollar | Archaeological prospection and remote sensing[END_REF]. These limitations confirm the potential usefulness of small, frequency-domain electromagnetic (also referred to as Slingram, looploop, or dipole-dipole) instruments, which can simultaneously measure the ground's apparent electrical conductivity and magnetic susceptibility [START_REF] Parchas | Simultaneous measurement of electrical conductivity and magnetic susceptibility of ground in electromagnetic prospecting[END_REF], and can now be fitted with a multi-receiver capability [START_REF] Saey | Electrical conductivity depth modelling with a multi-receiver EMI sensor for prospecting archaeological features[END_REF][START_REF] Bonsall | Assessment of the CMD Mini-Explorer, a new low-frequency multi-coil electromagnetic device, for archaeological investigations[END_REF].

The present paper focuses on the interpretation of measurements of the ground's in-phase magnetic susceptibility. Other measurements, such as that of the ground's magnetic viscosity [START_REF] Thiesson | TDEM magnetic viscosity prospecting using Slingram coil configuration[END_REF], or measurements combining both EMI and the magnetic technique, are reported to be highly advantageous in certain contexts [START_REF] Benech | Joint interpretation of E. M. and magnetic data for near-surface studies[END_REF][START_REF] Pétronille | Magnetic signal prospecting using multiparameter measurements : the case study of the Gallic site of Levroux[END_REF]. The aim of the present study is to assess the performance of a fast 3D linear interpretation algorithm, which is easy to implement on a laptop computer.

Through the use of the Moment Method (MoM), fast 3D inversion processes have already been proposed for DC resistivity prospection [START_REF] Brinon | Rapid 1D/3D inversion of shallow resistivity multipole data: examples in archaeological prospection[END_REF], and for the simultaneous interpretation of electrical conductivity and magnetic susceptibility EMI data [START_REF] Benech | Interpretation of shallow EMI resistivity and magnetic susceptibility measurements using rapid 1D/3D inversion[END_REF], whilst this processes is limited to reduced areas in the vicinity of targeted features. In the following, an inversion technique is proposed for the processing of magnetic susceptibility variation data derived from the MoM. Since the relative magnetic permeability of the soil varies only slightly (between 1.0001 and 1.01), it has been verified [START_REF] Tabbagh | The response of a tree dimensional magnetic and conductive body in shallow depth E.M. prospecting[END_REF]) that the Born approximation can be applied to the analysis of such data, thus paving the way for the use of linear inversion techniques: this is a sound approximation, and is systematically used in magnetic prospection when the so-called demagnetizing field is neglected [START_REF] Grant | Interpretation theory in applied geophysics[END_REF]West 1965, Scollar et al. 1990). On the other hand, it cannot be applied in the case of DC resistivity [START_REF] Dabas | 3-D inversion in subsurface electrical surveying-I: Theory[END_REF][START_REF] Buvat | Approximate 3D resistivity modelling using Fourier analysis of layer resistivity in shallow soil studies[END_REF] or conductivity electromagnetic (EM) methods because the soil's electrical conductivity varies over a wide range, between ≤ 0.1 and ≥ 100 mS/m. The implementation of a linear method has the advantage of being more efficient computationally, thus allowing the entire surveyed surface to be processed in one go.

Our proposed interpretation process is divided into two phases. In the first of these, the usual 1D non-linear inversion scheme is applied, point by point, to determine the vertical distributions of both electrical conductivity and magnetic susceptibility [START_REF] Zhang | Simultaneous reconstruction of 1-D susceptibility and conductivity from electromagnetic data[END_REF]. These vertical distributions are often sufficient to allow the terrain's structure to be characterized. If marked lateral changes require a 3D inversion the statistical mode of these distributions as a function of depth is taken to be the reference vertical profile.

Then, during the second phase, the interpreter determines the lateral variations in susceptibility with respect to this profile. Since the response of a contrasting magnetic body can be assumed to be linear, this process involves successively applying a 2D (horizontal plane) Fourier transform, non-linearly inverting the resulting spectrum of the 1D vertical distribution of magnetic susceptibility contrasting with respect to the reference profile, and finally computing the lateral susceptibility contrast variations by means of an inverse Fourier transform. The full process is illustrated in Figure 1.

Reminder of the 1D inversion process

The EMI instruments under consideration have a transmitting coil and at least one receiving coil at a metric separation expressed by L. The coils are located at a height d, above ground level. The coil orientations can be adjusted, and various possible configurations have been defined [START_REF] Frischknecht | Profiling methods using small sources[END_REF] with respect to the plane of the coils. The most commonly used configurations are HCP (horizontal coplanar), VCP (vertical coplanar), PARA (parallel -inclined at 55° from the horizontal plane), and PERP (perpendicular). The application of the instrument's primary field over a layered ground generates complex secondary fields, H s , which are expressed by Hankel Transform integrals. The full development of these expressions, and the definition of the subsequent approximations are described by [START_REF] Thiesson | Calibration of frequency domain electromagnetic devices used in near surface surveying[END_REF]. The method used to compute these transforms (i.e. the forward problem) is well known, and is similar to that used for vertical electrical sounding [START_REF] Ghosh | The application of linear filter theory to the direct interpretation of geoelectrical resistivity sounding measurements[END_REF]. The calculation results can be expressed in terms of apparent properties, this allows an initial approximate evaluation to be made of the soil's properties, thus simplifying comparisons between different instruments. The apparent magnetic susceptibility, κ a , is the susceptibility of a homogeneous terrain that would produce the same results with the same instrumental geometry (L, coil orientation, height above the surface).

The inverse problem is usually solved by starting from an a priori guess at the values of the unknown parameters, which are iteratively modified until a good fit between the computed results and the apparent experimental properties is reached. The number of parameters is limited by the number of different geometrical configurations which can be adopted by the instrument (coil orientations and L distances), since at the studied frequencies (VLF and LF) and depth ranges, the investigation depth is governed by the instrument geometry. Although a single parameter inversion can sometimes be useful (Guerin et al. 1996), several instruments (or a multi-coil configuration in the same apparatus) clearly provide a more detailed description of the terrain's structure.

It is important to note that with EMI data inversion, the vertical conductivity profile must first be determined (using the quadrature out of phase component of the responses), prior to inversion of the magnetic susceptibility profile because: (1) the conductivity distribution modifies the total magnetic field distribution in the ground, whereas the terrain's vertical susceptibility distribution has a negligible influence on the total magnetic field distribution [START_REF] Tabbagh | The response of a tree dimensional magnetic and conductive body in shallow depth E.M. prospecting[END_REF] and (2) when high, the conductivity may generate an in-phase response [START_REF] Thiesson | Calibration of frequency domain electromagnetic devices used in near surface surveying[END_REF]) that must be subtracted from the total in-phase response.

At the end of the 1D inversion step, the statistical mode of each model parameter is adopted, in order to define the reference vertical profile. A susceptibility response (d) κ a0 is associated with this profile, and is subtracted from the experimental response at each point (x, y) of the surveyed area, and for each receiver (at elevation d), thereby defining a contrasting magnetic susceptibility response d) y, (x, δκ a .

2D (x,y) inversion using a Fourier transform

Principle

Given δκ(x',y',z'), the contrasting magnetic susceptibility distribution inside the ground, and since the problem to be solved is linear, the response d) y, (x, δκ a that is added to the response of the reference profile can be expressed as a convolution product:

∫∫∫ ∞ - - - dz' dy' )dx' z' d y y x )IR(x z' y δκ(x = d) y, (x, δκ a , ' , ' , ' , ' (1)
,

where IR is the 'impulse response', the 3D Green function corresponding to the dipole source response to the instrumental configuration under consideration (analytical expressions corresponding to these functions are presented in more detail in [START_REF] Tabbagh | The response of a tree dimensional magnetic and conductive body in shallow depth E.M. prospecting[END_REF]). If the ground is discretized in the form of N successive layers having thicknesses e i centered at z i , and a magnetic susceptibility contrast δ

) z y κ(x i , ' , '
, equation ( 1) can be approximated by:

' ' ) , ' , ' ( ) , ' , ' ( ) , , ( 1 dy dx z d y y x x IR z y x e d y x i N i i i a - - - = ∑ ∫∫ = ∞ δκ δκ (2).
By applying a 2D Fourier transform to Eq. ( 2) one has:

) , , ( ) , , ( ) , , ( ˆ1 ∑ = - = N i i i i a z d v u I R z v u e d v u ) ) κ δ κ δ (3),
where a 'circumflex accent' indicates Fourier transformed functions, and (u, v) are the spatial frequencies corresponding to (x, y). Consequently, at each point (u, v) the expression for Eq. ( 3) corresponding to each level z i can be thought of as a single component of an Nequation linear system. If the number of different geometric coil configurations K is smaller than N, the system has no solution, however when K = N it can be solved directly. When K > N the system can be solved using the least squares method. Knowing δ ) z v, κ(u, i , it is straightforward to derive the solution for the problem of 3D susceptibility contrast via the inverse transformation of δ ) z y, (x, κ i ˆ.

The thin layer solution

Although commercial multi-configuration instruments are now commonly available and used in field surveys, a considerable volume of data has been acquired over the last forty years using just one type of configuration. It is thus important to reconsider and enhance the interpretation of this data. When a single configuration is used, only one layer can be considered and the resulting thin-layer interpretation is similar to that obtained with potential methods [START_REF] Grant | Interpretation theory in applied geophysics[END_REF]. In the case of magnetic prospection, this is a valuable solution [START_REF] Desvignes | About the determination of magnetic anomaly sources[END_REF]) since it allows: (i) the general pattern of the source body to be restored by correcting for the influence of the Earth's magnetic inclination on the shape of the anomaly; (ii) the maximum source depth to be assessed, and (iii) the magnetization contrast to be determined. Similar issues arise when interpreting EMI apparent magnetic susceptibility data: deformations introduced by the coil configuration, poor accuracy of depth assessments, and the contrast of any source body. This similarity has made it possible to use transformations between apparent susceptibility measurements and induced magnetization anomalies of the Earth's magnetic field, and has opened up significant perspectives for the simultaneous interpretation of both types of data [START_REF] Benech | Joint interpretation of E. M. and magnetic data for near-surface studies[END_REF]. In the present study, our analysis is limited to that of EMI measurements, for which the determination of the depth of the thin layer raises a specific issue: a relationship clearly does exist between the coil separation and the depth at which a susceptibility variation can be detected, if the sensitivity at a depth is too low, then the inversion near that depth will be poor.

We thus consider a series of synthetic cases, in which a 2m x 2m slab with a thickness of 0.2 m and susceptibility contrast equal to 200 10 -5 SI is displaced along the vertical axis.

The apparent susceptibility maps are inverted for each different position (change in depth) of the slab, and the quality of the inversion is assessed by comparing the susceptibility contrast produced by the inversion with the original value. Table 1(a) presents the results obtained for the different depths of the slab and coil separations using VCP configuration. Table 1(b) provides the results with the PERP configuration, Table 1(c) the results with the HCP configuration, Table 1(d) the results with the PARA configuration. Four main conclusions can be drawn for these comparisons:

(1) With the VCP configuration, the slab is correctly restored, even for a coil separation of L=2m and shallow slab depths, as well as for smaller coil separations and greater depths.

(2) With the PERP configuration, the results are comparable to those for the VCP configuration, and (with the exception of shallow slabs with VCP) the inverted values are generally 10% to 20% higher than the original values. For shallower slabs the outcome can probably be explained by the more complex lateral variations of the PERP dipole impulse response, and thus more problematic inversion of PERP data.

(3) With the HCP and PARA configurations the sign of the contrast is correctly restored (the apparent susceptibility anomaly is negative for greater d/L ratios with HCP) but not its magnitude (most often significantly amplified) (4) This preceding example, established using synthetic data, shows that good depths of investigation can be reached, even for the case of smaller values of L.

However, even for one thin layer, the simultaneous use of several configurations leads to slightly better results than the use of a single configuration.

Multi-layer inversion

At each point of the spectral domain, the solution δ

) z v, κ(u, i is represented by an N component vector ) κ , , (κ = κ δ N ... 1 r
, the solution for the linear equation system. While in practice limited to few components, the solution for the N linear equations system is confronted by an instability arising from the fact that an abnormal value in one layer can be compensated by another abnormal value, of the opposite sign, in another layer. The solution can be stabilized by adding an external constraint, for example by minimizing the norm of the solution vector. Although it is also possible to use this constraint to introduce the depth sensitivity dependence associated with each coil configuration, this sensitivity tends to be of the 'all or nothing' type, as shown in Table 1. Following several tests, two options were found to be potentially useful: either minimizing the square of the norm of the κ δ r vector (constraint

I): ∑ = N i Q 1 2 κ ) , (4),
or minimizing the sum of the squares of the differences between consecutive components (constraint II):

∑ - - = N i i Q 2 2 1 ) ( κ κ ) )
(5).

Thus, for K different instruments the minimized quantity is:

Q K j j theo a j a λ κ δ κ δ + - ∑ =1 2 , , exp, , ) ( ) ) (6),
where λ is chosen by the interpreter: for the second option the value λ=1 would be set, whereas with the first option it is more convenient to use the trace of the Jacobian matrix:

2 / trace = λ
. With constraint II the stabilization is stronger, the vertical variation of the magnetic susceptibility contrast minimum.

The application of this inversion procedure using synthetic data computed for PERP, VCP and PARA configurations and inter-coil distances varying between 0.7 and 1.5 m is presented in Table 2 for the case of constraint II and for λ=1. The body is 2m x 2m x 0.6m sized, it presents a 200 10 -5 SI susceptibility contrast, it is centered at z=0.5 m and divided in 5 layers (0.12m thick). Equivalently good results were obtained with reduced number of layers equal to: 3 (0.2 m thickness) or 4 (0.15 m thickness).

Tests with multi-body synthetic data

We consider a rectangular 25 m x 17 m area, meshed onto a regular 0.5m x 0.5m grid, containing four different magnetic bodies (Fig. 2) imbedded in a terrain with a homogeneous resistivity equal to 70 Ωm, and susceptibility equal to 30 10 -5 SI. The first body (A) is an Lshaped ditch with a 1m x 0.4 m cross section, centered at a depth of 0.4 m, the two perpendicular branches of which are 8 m in length and have a susceptibility contrast of 120 10 -5 SI. The second body (B) is a medium-sized square body with dimensions 2m x 2m x 0.75 m, centered at a depth of 0.6 m and having a susceptibility contrast equal to 100 10 -5 SI.

The third body (C) is a small superficial feature with dimensions 0.4m x 0.4m x 0.2 m, centered at a depth of 0.3 m depth and having a susceptibility contrast equal to 100 10 -5 SI.

The fourth body (D) is a slim feature with dimensions 0.4m x 0.4m x 1.6 m, centered at a depth of 1 m and having a susceptibility contrast equal to 100 10 -5 SI. The apparent magnetic susceptibility maps computed using MoM are shown in Fig. 3, for the seven different apparatus geometries described in Table 3. For all configurations, the transmitter-receiver line lies parallel to the x axis. The four features, with their different geometries, can be recognized in each map. These are characterized by ringing, especially in the case of the longest coil separations. It should be noted that the shortest apparatus provides the best description of the four magnetic bodies, in terms of the shape and the magnitude of the associated anomalies.

The seven datasets were inverted together, whilst considering three layers with different values of magnetic susceptibility, the first centered at z 1 =0.3 m with a thickness e 1 =0.2 m, the second centered at z 2 =0.5 m with a thickness e 2 =0.2 m, and the third centered at z 3 =0.8 m with a thickness e 3 =0.4 m. When 1D inversion is applied (Fig. 4), although the four features are correctly identified, this inversion fails to correct for their deformed shapes and for the apparent anisotropy, produced by the different coil orientations. This issue is particularly pronounced for small features, whose precise locations are difficult to restore. In addition, ringing with sign changes cannot be eliminated, and is even amplified in the third layer in which there should be no variations corresponding to the superficial features (A) and (C). These spurious variations are a consequence of the instrument's low sensitivity to 3D changes at the depth of these features.

The results produced by 3D linear inversion are presented in Fig. 5, showing that the body's shapes are faithfully restored, and the susceptibility values are in good agreement with the original values for the first two layers. However, in the case of feature (A) the third layer is again characterized by greater amplitude variations, resulting from the fact that the instruments' sensitivities are too low to constrain the solution.

Interpretation of field data

In order to assess the potential of this approach for data inversion under field conditions, we analyzed the data acquired using various instrument geometries, at four different archaeological sites characterized by different climatic conditions and soil environments. In all cases we started from the apparent susceptibility data sets which have been obtained after checking the instrument calibration [START_REF] Thiesson | Calibration of frequency domain electromagnetic devices used in near surface surveying[END_REF]) and transforming the in-phase secondary field measurements. Where the electrical conductivity is high (Medamud case), the apparent electrical conductivity was first determined using the quadrature response and then used to calculate the in-phase response which was algebraically subtracted from the total in-phase response.

Neolithic enclosure at Balloy (Seine et Marne, France)

The eastern section of this middle-neolithic 'Passy' type of funeral enclosure [START_REF] Mordant | Le complexe des Réaudins à Balloy: enceinte et nécropole monumentale. La Culture de Cerny[END_REF]) has been the object of multi-method tests. This enclosure was detected using both electrical (1 m square array) and SH3 prospection (see appendix 1 for the characteristics of the latter device), but was not detected by magnetic prospection using a fluxgate gradiometer with 1 nTm -1 sensitivity [START_REF] Hesse | Rapport de synthèse sur les prospections géophysiques des vestiges archéologiques: Direction régionale des antiquités d'Ile de France[END_REF]. This failure was explained by the use of a 3D nonlinear inversion (Bénech et al. 2016), applied to a selected zone (indicated by a rectangular outline in Fig. 6), showing that the SH3 measurements had revealed a 1.4m x 0.4 m crosssectional feature, surrounded by gravel, with magnetic susceptibilities of respectively 51 x 10 - 5 and 20 x 10 -5 SI. Using these parameters, the induced magnetization anomaly determined for a fluxgate vertical gradiometer is less than 0.5 nTm -1 , and even with the addition of viscous magnetic remanent magnetization the anomaly cannot clearly overpass 1nTm -1 .

In the present study we consider just one layer (since only one instrumental dataset was recorded), of 0.4m thickness and centred at a depth of 0.45m. For the same selected zone, the following results are obtained:

(1) as could be expected, the 1D inversion produces the same image as that generated using the apparent susceptibility measurements (Fig. 6b),

(2) the susceptibility variations are in agreement with the 3D non-linear inversion, with a magnetic susceptibility of 40 10 -5 SI for the ditch filling, in contrast with a value of 10 x 10 -5 SI for the immediately surrounding gravel. This calculation was made for the totality of the 26m x 26 m surface shown in Fig. 6, and required just 7s of computing time, whilst using the same laptop computer (4 Go RAM, 2.5 GHz) the non-linear inversion over the 2.5m

x 5m area required 267s of computing time.

Neolithic settlement of Perdika 2 (Central Greece)

This neolithic settlement, in this region referred to as a magoula, is located on the eastern Thessalian plain of Greece and was surveyed using different methods in the frame of the IGEAN project (Innovative Geophysical Approach for the study of Early Agriculture villages of Neolithic Thessaly) [START_REF] Kalayci | How efficient is an integrative approach in archaeological geophysics, Comparative case studies from neolithic settlements in Thessaly (central Greece)[END_REF]. The magnetic survey revealed a complex system of enclosures (Fig. 7). In order to acquire a better description of what is thought to be an enclosure entrance, a 40m x 40m area was prospected using the CMD instrument (see Appendix I for its characteristics). The ditches are more clearly discernible by their magnetic susceptibility contrast than by their electrical conductivity. The resistivity variations are relatively small, lying within a moderate range of values: a 150 Ωm layer of topsoil, and subsoil variations ranging between 30 and 100 Ωm. As the magnetic susceptibility data acquired with the shorter coil separation (L=0.32m) was too noisy for suitable interpretation, we made use of four data sets for the inversion (Fig. 8): HCP with 0.71 m and 1.18 m separations, and VCP with 0.71 m and 1.18 m separations. The topsoil layer was assumed to be homogeneous. In order to assess the vertical extent of the archaeological features, we chose to divide the subsoil into two contiguous layers of 0.5 m thickness, with the first lying between 0.2 and 0.7 m, and the second between 0.7 and 1.2 m. The images produced by point-by-point 1D inversion and by 3D linear inversion of the data are presented in Figs. 9a and9b, respectively. In the 1D inversion, the second layer reveals some features that are different to those found in the first layer, which is characterised by strong susceptibility variations. Also performed using constraint I, see equation ( 4), the 3D linear inversion shows that the ditches are shallow, and that different features appear at greater depths, but with smaller variations in susceptibility.

Middle Kingdom Kôm of Medamud (Egypt)

The In the present study we consider a 44m x 60m zone from this site. Fig. 10 shows three apparent susceptibility maps, produced using inter-coil distances of 0.32, 0.71 and 1.18 m, respectively. Fig. 11 shows the results of the inversions when two layers are considered, with the first situated between 0.05 and 0.35 cm, and the second between 0.35 and 0.65 m: Fig. 11a and Fig. 11b show the results of the 1D and 3D inversions (both using constraint I), respectively. Both inversions show that the magnetic features have a significant vertical extent (which probably exceeds the range of the instrument), and that their susceptibilities are high. However, comparisons between these different results underscore the fact that when 1D inversion is used, the second layer has higher values, in the range [350 -1250] 10 -5 SI, whereas the first layer remains close to a modal value of 200 10 -5 SI, which is the same as that found for both layers with the 3D inversion. Again in this case, the second layer of the 3D inversion is characterised by a narrower range of variations ([180 -225] 10 -5 SI) than the first layer ([130 -330] 10 -5 SI).

These maps clearly show the presence of a mud brick wall which corresponds to lower apparent susceptibilities. By comparing these results with measurements (MS2D Bartington Ltd) taken over another feature we obtain values in the range [60 -90] 10 -5 SI for the mud bricks and a modal value of 178 10 -5 SI for the surrounding magnetic soil (Relats-Montserrat et al. 2016(b)). When applying a full 3D inversion to the rectangular 12m x 3 m area (marked in Figure 11b) we obtain, by considering a wall of 0.6m thickness and 4m width, a susceptibility contrast of -92 10 -5 SI with the surrounding soil. These results are in fair agreement with the results mapped in Figure 10b and with the 3D linear inversion using one layer of 0.6m thickness which gives a -100 10 -5 SI susceptibility contrast.

Destroyed medieval city of Thérouanne (Pas-de-Calais, France)

Thérouanne was originally a Gallic, then an important Gallo-roman city, and during medieval times was one of the main centres of the Picardie region in France. In 1553 it was besieged and totally levelled by Charles Quint. This location was settled again as a small village, only at the end of the XIX th century. The aim of the survey was to identify the remains of the ancient medieval city, through a series of EMI and electrical prospection campaigns. The test area considered here covers a surface area of 180m x 40m, and was surveyed using the DualEM instrument (see Appendix I) in both HCP and VCP base configurations. From this archaeological survey, we produced six different apparent magnetic susceptibility maps, derived from the instrument's in-phase recordings made with the HCP 1m, HCP 2m, PERP 1.1m, PERP 2.1m, VCP 1m and VCP 2m configurations. In the present analysis, the longest inter-coil separations of 4 m and 4.1m were not used, and the HCP 1m

and PERP 2m channels produced corrupted data. Thus, only four independent magnetic susceptibility maps could be used, as shown in Fig. 12.

The electrical resistivity maps of the terrain reveal a 70 Ω m topsoil layer, followed by a conductive subsoil layer with a resistivity between 10 and 60 Ωm, suggesting a significant clay content (the salt spread under the order of Charles Quint was probably leached out well before the time of the prospection campaign). The apparent magnetic susceptibility is generally high, and higher for the PERP 1m and HCP 2m measurements than for both VCP configurations (different susceptibility scales are used in Fig. 12). The change in sign of the susceptibility response for the HCP 2m measurement, at approximately 0.75m, suggests that most of the observed features have shallow locations. Three layers were considered for the inversions: the first between 0.25 and 0.75 m; the second between 0.75 and 1.25 m; and the third between 1.25 and 2.25 m. The results of the 1D point-by-point inversion are presented in Fig. 13, and those of the 3D linear inversion are shown in Fig. 14. Constraint I was applied for both inversions. Although the mean values of susceptibility decrease slightly with depth in the 1D inversion, the lateral variations have the same aspect and reproduce those of the four apparent susceptibility maps. A large, 40 m diameter semi-circular feature, centred 18m to the north of the image centre, can be seen in all three layers, thus suggesting that it has a substantial vertical extent. In the 3D inversion, the first layer is clearly more magnetic than the remaining two layers (the scales are different in Fig. 14), and the large semi-circular feature is less noticeable in the deepest layer.

Conclusion

The new linear inversion technique presented here has been tested on synthetic and field data. When compared with the original model, or with the complete 3D MoM inversion technique, it provides reliable results for the shape of the sought features, as well as the magnitude of the contrasts. It is shown to be highly efficient for the correction of anomalous deformations, caused by the instrument's geometry, of a feature's geometrical outlines. By using the instrument's dipole impulse response in the analysis, this technique takes the influence of the real instrumental geometry into account, including that of the coil orientations. Its results are thus significantly different to those obtained with 1D point-topoint inversions, in terms of susceptibility variations as a function of depth, making it possible to improve the identification of the vertical boundaries of a given feature.

It should be underlined, however, that this technique is limited by the number of different available instrument geometries: in practice the magnetic susceptibility contrasts of only a rather small number of layers at a limited number of depths can be inversed. It is very valuable to have even this limited information on depth and susceptibility and where necessary the prospector has the possibility to increase the number of instruments geometries by considering measurements acquired at several altitudes with the same device.

Compared to 3D complete MoM inversion, the rapidity of this technique not only leads to a significant gain in time, but also provides the interpreter with the ability to test several stabilization process options, such as the choice of the number of layers or layer thicknesses. This also makes it possible to assess the probable vertical extent of the different features, through the use of a vertically translated thin layer. Figure 5: Results of the 3D linear inversion when three layers are considered (using the same geometries as in Fig. 3), and the norm of the κ vector is minimized using λ=trace/2

Figure captions

  site of Medamud is located 4 km to the north of Karnak, and was occupied for a long period of time, from the Middle Kingdom of Egypt until the Roman era (Relats-Montserrat 2016(a)). The geophysical survey initiated by the IFAO tended to characterize a Middle Kingdom pottery workshop area. The survey was carried out with the CMD (see Appendix I) instrument in a VCP configuration, in order to simultaneously describe the soil's electrical conductivity and magnetic susceptibility. In the area under consideration the finegrained, thick archaeological layer present just below the surface has a very low resistivity, centred on 16 Ωm, for which a correction is required to eliminate the in-phase component of the conductivity response

Figure 1 :

 1 Figure 1: Flowchart demonstrating the full process with both a 1D non-linear inversion and a

Figure 2 :Figure 3 :

 23 Figure 2: Description of the different bodies used to generate synthetic data

Figure 4 :

 4 Figure 4: Results of the 1D inversion when three layers are considered: the first is centered at

Figure 6 :

 6 Figure 6: Balloy neolithic enclosure: (a) SH3 apparent magnetic susceptibility, the relatively high level is explained by the high susceptibility of the upper layer of soil (approx. 100 x 10 -5

Figure 7 :

 7 Figure 7: Perdika2 (Central Thessaly, Greece) magnetic map: pseudo-gradient of the vertical

Figure 8 :

 8 Figure 8: Perdika2 site, apparent magnetic susceptibility maps measured using the CMD

Figure 9 :

 9 Figure 9: Perdika2 site, (a) inverted magnetic susceptibility of the 0.2 -0.7 m and 0.7 -1.2 m

Figure 10 :

 10 Figure 10: Medamud site (Egypt), three apparent magnetic susceptibility maps measured

Figure 11 :

 11 Figure 11: Medamud site, (a) inverted magnetic susceptibility of the 0.05 -0.35 m and 0.35 -

Figure 12 :

 12 Figure 12: Thérouanne (Pas-de-Calais, France) site, four apparent magnetic susceptibility

Figure 13 :

 13 Figure 13: Thérouanne site, inverted magnetic susceptibility of the 0.25 -0.75 m, 0.75 -

Figure 14 :

 14 Figure 14: Thérouanne site, inverted magnetic susceptibility of the 0.25 -0.75 m, 0.75 -

  Fig.1

  Fig. 2

Fig. 3

 3 Fig.3

Fig. 4

 4 Fig. 4

  Fig.5

  Fig.6

Fig. 9a

  Fig. 9a

  Fig. 9b

Fig. 11a

  Fig. 11a

Fig. 11b

  Fig. 11b

Fig. 12

 12 Fig. 12

  Fig.13

Fig. 14

 14 Fig. 14

Table captions Table 1 :

 captions1 Inverted susceptibility contrast (in 10 -5 SI) as a function of the slab's depth of the

	slab, z, and the separation between transmitting and receiving coils, L (the original value
	being 200 10 -5 SI), (a) VCP configuration, (b) PERP configuration (the symbol X means that
	the shape of the slab is not correctly restored and unidentifiable), (c) HCP configuration and
	(d) PARA configuration

Table 2 :

 2 Inverted susceptibility contrast when a 2m x 2m x 0.6 m body, of 200 10 -5 Si susceptibility contrast is divided into 5 thin layers of 0.12m thickness, determined using respectively five, six and seven different instrumental configurations (constraint II is used in the inversion).

Table 3 :

 3 Characteristics of seven instrumental configurations used to interpret inversionprocesses with synthetic data.

		L=0.5 m	L=0.7 m	L=1 m	L=1.5 m	L=2 m
	z=0.3 m	195	187	153	169	153
	z=0.5 m	227	224	209	226	200
	z=0.7 m	245	235	221	236	234
	z=1 m	242	241	239	241	229
	z=1.5 m	250	226	219	249	205
	Table 1(a)
		L=0.5 m	L=0.7 m	L=1 m	L=1.5 m	L=2 m
	z=0.3 m	206	258	183	X	X
	z=0.5 m	235	244	217	259	219
	z=0.7 m	242	240	233	253	228
	z=1 m	252	245	252	248I	230
	z=1.5 m	X	X	X	249	193
	Table 1(b)
		L=0.5 m	L=0.7 m	L=1 m	L=1.5 m	L=2 m
	z=0.3 m	296	120	240	225	X
	z=0.5 m	271	221	215	233	302
	z=0.7 m	266	304	257	224	284
	z=1 m	271	282	263	292	276
	z=1.5 m	277	279	249	280	248
	Table 1(c)

Appendix I Characteristics of the devices used for field studies SH3: this apparatus is a laboratory prototype designed in 1977 [START_REF] Parchas | Simultaneous measurement of electrical conductivity and magnetic susceptibility of ground in electromagnetic prospecting[END_REF].

It has a PARA coil orientation (the two coils have parallel axes aligned at 35° from the vertical, such that their direct coupling is null in free space), a 1.5 m coil separation, a coil center height of d=0.2 m, and is operated at 8.04 kHz.

CMD: this apparatus is a multi-receiver EMI (Gf Instruments, Ltd, Brno) comprising one transmitter coil and three receiver coils, located at distances of 0.32, 0.71 and 1.18 m from the transmitter. All the coils are coplanar, allowing the instrument to be used in either HCP or VCP configurations. The instrument's operating frequency is 30 kHz, and it can be used in a continuous recording mode by a mobile operator in the field, with an above-ground clearance of d=0.12 m. DualEM 421S: this apparatus is a multi-receiver EMI (DualEM Ltd, Milton) operated at 9 kHz. It associates one horizontal transmitter loop with three pairs of receivers. In each pair, the first receiver is horizontal, allowing HCP measurements to be made. By rotating the entire apparatus, VCP configuration measurements can also be made. The second receiver of each pair is oriented in a radial direction from the transmitter, allowing PERP configuration measurements to be used. The receivers of the first pair are located at respectively 1m and 1.1m from the transmitter, those of the second pair at 2m and 2.1 m, and those of the third pair at 4m and 4.1m. However, in the archaeological surveys presented here, data from the third pair were not considered. Data from the VCP 1m, VCP 2m, HCP 1m, HCP 2m, PERP 1.1m and PERP 2.1m are used in the present study.


Catherine Gérikas Ribeiro
email: catherine.gerikas@gmail.com

Dominique Marie

Adriana Lopes

Frederico Santos

Daniel Pereira Brandini

Vaulot

Catherine Gérikas Ribeiro

Adriana Lopes Dos Santos

Frederico Pereira Brandini

Daniel Vaulot

Estimating microbial populations by flow cytometry: Comparison between instruments

Keywords: flow cytometry, microbial oceanography, heterotrophic bacteria, picoplankton, nanoplankton

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Flow cytometry (FCM) is a well-established technique [START_REF] Marie | Enumeration and cell cycle analysis of natural populations of marine picoplankton by flow cytometry using the nucleic acid stain SYBR Green I[END_REF][START_REF] Gasol | Using flow cytometry for counting natural planktonic bacteria and understanding the structure of planktonic bacterial communities[END_REF] used since the 1980s [START_REF] Trask | Analysis of phytoplankton by flow cytometry[END_REF][START_REF] Olson | Marine phytoplankton distributions measured using shipboard flow cytometry[END_REF] for enumeration and characterization of marine micro-organisms. FCM analysis of planktonic communities fulfills the scientific demands of rapid and accurate cell counting, as it considerably reduces the bias introduced by visual counting [START_REF] Marie | Phytoplankton Cell Counting by Flow Cytometry[END_REF]. By simultaneously recording several parameters during analysis, FCM allows the discrimination of pico-and nanoplankton populations and the estimation of their abundance, cell size, and pigment content [START_REF] Marie | Phytoplankton Cell Counting by Flow Cytometry[END_REF], both by natural (chlorophyll, phycoerythrin) or induced (fluorescent dyes) fluorescence [START_REF] Marie | Enumeration and cell cycle analysis of natural populations of marine picoplankton by flow cytometry using the nucleic acid stain SYBR Green I[END_REF].

The flow cytometer registers events as cells are aligned in a fluid stream and flow through a beam of focused light usually provided by one or several lasers. For each particle, scattered light and emitted fluorescence are converted to digital signals and recorded. A flow cytometer comprises three main systems: fluidics (particle transport), optics (laser beam and optical filters), and electronics (signal conversion into electronic data). Detectors for scattered light located at 180° and 90° from the light source are called forward scatter (FSC) and side scatter (SSC), respectively. Fluorescence at different wavelengths (typically green, orange and red) is also recorded. Signals associated with each parameter are displayed as cytograms, which are used to discriminate and count different populations based on scattering and fluorescence features. Phytoplankton populations can be differentiated by FCM according to specific values of the recorded parameters (FSC, SSC, red or orange fluorescence).

Bacteria are in general detected after staining with a nucleic acid strain such as SYBR Green-I [START_REF] Marie | Enumeration and cell cycle analysis of natural populations of marine picoplankton by flow cytometry using the nucleic acid stain SYBR Green I[END_REF]. Two different groups can be distinguished based on their apparent nucleic acid content (differences in fluorescence intensity) and side scatter signal (SSC): high nucleic acid (HNA) and low nucleic acid (LNA) bacteria. The function and ecological importance of these two groups is far from being fully understood [START_REF] Bouvier | A comparative study of the cytometric characteristics of High and Low nucleic-acid bacterioplankton cells from different aquatic ecosystems[END_REF][START_REF] Van Wambeke | Vertical and longitudinal gradients in HNA-LNA cell abundances and cytometric characteristics in the Mediterranean Sea[END_REF], and several studies have addressed these nucleic acid content differences in terms of ecological traits, such as bacterial activity and production [START_REF] Morán | Assessing the relevance of nucleic acid content as an indicator of marine bacterial activity[END_REF][START_REF] Ortega-Retuerta | Exploring the relationship between active bacterioplankton and phytoplankton in the Southern Ocean[END_REF][START_REF] Van Wambeke | Vertical and longitudinal gradients in HNA-LNA cell abundances and cytometric characteristics in the Mediterranean Sea[END_REF].

Two main groups of autotrophic prokaryotes dominate picoplankton: Prochlorococcus and Synechococcus. Prochlorococcus is ubiquitous in the euphotic zone of tropical oceans, being considered the most abundant photosynthetic organisms on the planet (Partensky et al. 1999b), and its discovery was only made possible with the development of flow cytometry [START_REF] Chisholm | A novel free-living prochlorophyte abundant in the oceanic euphotic zone[END_REF].

Prochlorococcus is discriminated by its small scattering and low red fluorescence (chlorophyll).

Synechococcus is widely distributed in marine environments, being particularly abundant in well-lit and nutrient rich top layers of the oceans (Partensky et al. 1999a). One of the key parameters that allows Synechococcus populations to be discriminated by FCM is the phycoerythrin content (orange fluorescence). Different Synechococcus clades can show distinct fluorescence signatures [START_REF] Olson | Pigments, size, and distribution of Synechococcus in the North Atlantic and Pacific Oceans[END_REF][START_REF] Thompson | A multi-laser flow cytometry method to measure single cell and population-level relative fluorescence action spectra for the targeted study and isolation of phytoplankton in complex assemblages[END_REF], as a result of different phycobilisome composition [START_REF] Scanlan | Ecological genomics of marine picocyanobacteria[END_REF]). Pico-and nanoeukaryotes are important contributors to global primary productivity [START_REF] Li | Primary production of prochlorophytes, cyanobacteria, and eucaryotic ultraphytoplankton: Measurements from flow cytometric sorting[END_REF], and due to their larger cell size, they often contribute to an important share of autotrophic biomass in the oceans [START_REF] Zubkov | Picoplanktonic community structure on an Atlantic transect from 50°N to 50°S[END_REF]. Picoeukaryotes, which cells range from 0.8 µm to 2-3 µm [START_REF] Simon | Characterization of oceanic photosynthetic picoeukaryotes by flow-cytometry[END_REF], present well-defined cytometric signatures by FCM, while nano-eukaryotes populations are less well defined.

FCM analysis has led to numerous advances in marine microbial ecology, although cost and maintenance expenses were prohibitive for many laboratories until recently [START_REF] Gasol | Using flow cytometry for counting natural planktonic bacteria and understanding the structure of planktonic bacterial communities[END_REF][START_REF] Vives-Rego | Current and future applications of flow cytometry in aquatic microbiology[END_REF]. Since the first cytometry-based field study made by [START_REF] Olson | Marine phytoplankton distributions measured using shipboard flow cytometry[END_REF], on-board flow cytometry has become a crucial tool in the investigation of both autotrophic and heterotrophic picoplanktonic communities [START_REF] Legendre | Flow cytometry in oceanography 1989-1999: Environmental challenges and research trends[END_REF]. The manufacturing of low cost compact benchtop flow cytometers such as the BD Accuri™ C6, the Millipore Guava ® or the Applied Biosystems Attune ® has facilitated the use of FCM to study of phytoplankton communities around the world, due to easy handling, automatic sampling and easy transportation (a critical quality for field measurements). However, these low cost instruments can be less sensitive than laboratory based flow cytometers, due to less sophisticated optical and/or electronic systems.

A lower sensitivity is usually not a problem for bacteria which are detected after staining with strongly fluorescing dyes such as SYBR Green, or for small eukaryotes whose pigment content is relatively high. However, this is not the case for cyanobacteria like Prochlorococcus, for which the concentration of photosynthetic pigments per cell is as much as 50-100 fold lower in cells exposed to high light as a result of photo-acclimation [START_REF] Sosik | Chlorophyll fluorescence from single cells: Interpretation of flow cytometric signals[END_REF][START_REF] Olson | Pigments, size, and distribution of Synechococcus in the North Atlantic and Pacific Oceans[END_REF]), creating 'dim' populations in the surface layers. Such low fluorescence explains why Prochlorococcus escaped detection by researchers using epifluorescence microscopy or even during the first use of FCM on board oceanographic ships [START_REF] Olson | Marine phytoplankton distributions measured using shipboard flow cytometry[END_REF]. To overcome the problem of low sensitivity flow cytometers, both direct and indirect approaches to infer Prochlorococcus abundance have been developed, such as changes in cytometer optical set up to improve excitation energy or fluorescence detection [START_REF] Dusenberry | Increasing the sensitivity of a FACScan flow cytometer to study oceanic picoplankton[END_REF]Partensky et al. 1999b) and the use of mathematical corrections [START_REF] Zubkov | Picoplanktonic community structure on an Atlantic transect from 50°N to 50°S[END_REF][START_REF] Crosbie | Abundance , distribution and flow-cytometric characterization of picophytoprokaryote populations in central (17 o S) and southern (20 o S shelf waters of the Great Barrier Reef[END_REF].

In this paper, we compare data obtained on marine microbial communities with two flow cytometers, the FACSCanto™ and the Accuri™ C6 (hereafter named as CANTO and C6). Although manufactured by the same company (BD Biosciences, San Jose, CA, USA), these cytometers present distinct fluorescence excitation/detection technical features (Table 1). Differences in laser, optics and electronic systems can potentially affect sensitivity and resolution, influencing the accuracy of field measurements. We analyzed heterotrophic marine bacteria, photosynthetic eukaryotes and cyanobacteria on a set of marine samples from the South Atlantic Ocean (displaying both nutrient and light gradients within the water column). While both instruments produced equivalent data for bacteria and eukaryotes, cyanobacteria, especially Prochlorococcus, were severely under-estimated with the C6 instrument, and procedures previously suggested to correct the data proved ineffective.

Materials and methods

Sampling

Water samples were collected onboard the R/V "Alpha Crucis", between 31/10/2013 and 23/11/2013. The surveyed area was located between latitude 23°11'S -30°52'S and longitude 39°22"W -49°09"W, along 3 transects (TR1, TR2 and TR3), in the South West Atlantic off Brazil, reaching the 3510 meters isobath (Fig. 1). The sampling strategy comprised cross-shelf transects with 5 depths per station within the euphotic zone for TR1 and TR2, as well as 12 surface samplings for TR3, for a total of 102 samples. Three water masses were sampled during the cruise: the warm and oligotrophic Tropical Water, the cold and nutrient rich South Atlantic Central Water and the Coastal Water, with highly variable features [START_REF] Castro | Physical oceanography of the western Atlantic continental shelf located between 4 N and 34 S[END_REF]. Except for TR3 samples, which were collected using a polycarbonate bucket, seawater samples were collected with 12 L Niskin bottles attached to a Seabird ® CTD-rosette system (Sea-Bird Electronics, Bellevue, WA, USA), divided into cryotubes, preserved with 0.1% glutaraldehyde, final concentration (modified from [START_REF] Vaulot | A simple method to preserve oceanic phytoplankton for flow cytometric analyses[END_REF], incubated for 10 minutes in the dark, flash-frozen in liquid nitrogen and stored at -80°C until analysis.

Flow cytometry analysis

Two flow cytometers were used in this study: a BD FACSCanto II™ and a BD Accuri™ C6 (Table 1). Samples were counted simultaneously on both cytometers located in the same room, in order to avoid any possible bias by manipulation or time span between measurements. The tubing of the C6 was new and fluidics were calibrated for precise volume measurements as recommended by the manufacturer (Section 4.13 of manual).

Samples were first analyzed unstained to enumerate phototrophs. Fluorescent beads (0.95 G Fluoresbrite ® Polysciences, Warrington, PA) were added in each sample in order to normalize parameters [START_REF] Marie | Enumeration and cell cycle analysis of natural populations of marine picoplankton by flow cytometry using the nucleic acid stain SYBR Green I[END_REF]. A second analysis was performed to enumerate heterotrophic bacteria after staining with SYBR Green ® (1:10000, final concentration) (Ref-S7585, Life Technologies, Eugene, Oregon).

On the C6, for enumerating phytoplankton, 200 µL of sample were analyzed at the "high" rate (66 µL.min -1 ) with a threshold set at 700 on red fluorescence (FL3-H). To enumerate heterotrophs, 60 µL of SYBR Green stained samples were run at "medium" rate (35 µL.min -1 ) and the threshold was set at 700 on green fluorescence (FL1-H). In both cases, thresholds were determined by running 0.2 µm filtered sea water sample and lowering the values until electrical or optical noise appears.

On the CANTO for enumerating phytoplankton, samples were run for 3 min with a rate of 72 µL.min -1 and with the discriminator set on red fluorescence at 200. For bacterial enumeration, SYBR Green stained samples were run for 2 min at a rate of 60 µL.min -1 and the threshold was set on green fluorescence at 500. Flow rate was determined by the method described by [START_REF] Marie | Enumeration and cell cycle analysis of natural populations of marine picoplankton by flow cytometry using the nucleic acid stain SYBR Green I[END_REF]. A known volume of seawater was injected on the CANTO for at least 10 min. Then the remaining volume is measured and the rate is determined by dividing the difference between initial and final volumes by the injection time.

Data were analyzed with the Flowing Software ® 2.5 (http://www.flowingsoftware.com). Each population was identified on the cytograms on the basis of its scatter and fluorescence signals (Fig. S1).

Each parameter was normalized to that of the reference beads (0.95 µm). Cell counts for each red fluorescence value were exported from the single parameter histogram. The resulting spreadsheet (File S1) was used in subsequent analysis with the R software (R Development Core Team, 2013).

Near the surface the red chlorophyll fluorescence of the picophytoplankton decreases due to photoacclimation [START_REF] Partensky | Photoacclimation of Prochlorococcus sp. (Prochlorophyta) Strains isolated from the northern Atlantic and the Mediterranean sea[END_REF][START_REF] Dusenberry | Photoacclimation kinetics of single-cell fluorescence in laboratory and field populations of Prochlorococcus[END_REF][START_REF] Kulk | Distinct differences in photoacclimation potential between prokaryotic and eukaryotic oceanic phytoplankton[END_REF]. Therefore, for a fraction or even all of the Prochlorococcus and Synechococcus populations, fluorescence can fall below the detection threshold (Fig. 2). For the case where only a part of the population was in the noise, we modified the correction procedure described by [START_REF] Crosbie | Abundance , distribution and flow-cytometric characterization of picophytoprokaryote populations in central (17 o S) and southern (20 o S shelf waters of the Great Barrier Reef[END_REF] and implemented it as an R routine (File S2). This correction assumes that the red fluorescence distribution of these populations has a log-normal shape [START_REF] Crosbie | Abundance , distribution and flow-cytometric characterization of picophytoprokaryote populations in central (17 o S) and southern (20 o S shelf waters of the Great Barrier Reef[END_REF][START_REF] Shapiro | Practical flow cytometry, 3rd Edition[END_REF] and that, when the left part of the distribution is partially in the noise, the left part can be extrapolated from the right part. The R routine takes as input histograms produced by the Flowing Software ® (but can be adapted to other data formats) and outputs uncorrected and corrected cell abundance data (output data examples can be found in Files S2 andS3). Three cases can occur for a given population.

1. The mode of the histogram is not visible (e.g. Fig. 2G). In this case, the population is considered to be mostly within noise, without the possibility of counting or correction. Hence, this population is removed from the dataset and labeled as 'cells in noise' by the R routine.

2. The mode is visible, but the left part of the distribution is partly below the noise level (e.g. Fig. 2A). In this case, the abundance of each population is calculated as the double of the right part of the histogram, from the mode to the maximum, and samples are labeled as 'correction' 3. The mode is visible and the distribution is totally out of the noise (e.g. Fig. 2C). In this case no correction is performed and the initial output value is kept ('no correction' samples).

In some cases, the automatic correction needs some degree of visual confirmation, especially for deeper samples with low cell numbers resulting in noisy histograms (e.g. Fig. 2Q or 2W).

Therefore, the R routine provides a graphical output of the histogram for each sample (File S3), allowing the user to visually confirm whether the automatic labeling (case 1, 2, or 3 above) is correct.

Statistical analyses were performed with the PRISM ® 7 software (http://www.graphpad.com/scientificsoftware/prism).

Results

Bacterial populations were well resolved for both CANTO and C6 for all samples (Fig. S1).

High Nucleic Acid and Low Nucleic Acid bacterial populations were consistently distinguished with both cytometers and there was a very good correlation for both HNA and LNA (R 2 = 0.85 and 0.89, respectively) between the data obtained on the different instruments (Fig. 3A,B). For HNA slope was statistically different from 1 (p < 0.0001) and abundances estimated by the C6 were consistently lower by 10-15% compared to the CANTO.

The chlorophyll fluorescence of Prochlorococcus and Synechococcus decreases from the deeper layers to the surface in response to photoacclimation. For samples near the surface, cells from both populations can be partly or totally in the noise depending on the instrument sensitivity. Fig. 2 demonstrates clearly that the C6 is less sensitive than the CANTO by at least a factor of 10 (note for example Fig. 2C and 2I for Prochlorococcus at 110 m for the CANTO and C6, the mode of the histogram is in the noise for the C6 and about 10 times higher than the noise level for the CANTO).

With the C6, Prochlorococcus populations were completely in the noise above 100 m (55 out of 102 samples, Fig. 2 and4, Table 2) and for Synechococcus partly or completely in the noise above 50 m (15 out of 102 samples, Fig. 2 and4, Table 2). With the CANTO, only surface Prochlorococcus were partly in the noise and Synechococcus was always fully resolved (Fig. 2 and4, Table 2). When cells were only partly in the noise (i.e. when the histogram mode was clearly visible, e.g. Fig. 2U), we estimated the part of the population that was in the noise using the approach proposed by [START_REF] Crosbie | Abundance , distribution and flow-cytometric characterization of picophytoprokaryote populations in central (17 o S) and southern (20 o S shelf waters of the Great Barrier Reef[END_REF] (see Material and Methods section). The comparison between the C6 data which required correction and the CANTO data that did not require correction allowed us to assess the validity of this approach (Fig. 3C andD).

Clearly some data points that are corrected (grey squares) appear as outliers and are severely underestimated with the C6, even after correction (Fig. 3C andD). The slopes for non-corrected samples (solid circles) are significantly different from 1 (p < 0.0001), being respectively 0.75 (Prochlorococcus) and 0.77 (Synechococcus) (Fig. 3C andD), which corresponds to 25% lower abundance on average with the C6. Vertical cross sections of two transects (Fig. 4) illustrate that, while the CANTO provides fully resolved vertical profiles for both Prochlorococcus and Synechococcus population, the data from the C6 cannot be used in the upper layer (roughly from 100 m to the surface).

Pico-and Nano-eukaryotes were always above the detection limit for both instruments with excellent correlation between the two instruments (R² = 0.94 and 0.69, respectively, Fig. 3E andF). As for the other populations, the slopes were significantly different from 1 (p < 0.0001). While picoeukaryotes were about 15% more abundant with the C6, it was the reverse for nanoeukaryotes which were slightly underestimated by the C6.

Discussion

The analysis of planktonic communities by flow cytometry is complex because the distinctive cell features of each population may change with depth, diel cycle and nutrient conditions (e.g. Vaulot and Marie 1999). From the six planktonic groups studied here, four were well resolved by both the CANTO and C6 flow cytometers: HNA and LNA heterotrophic bacteria, autotrophic pico-eukaryotes and nano-eukaryotes. Abundance measured by both instruments were tightly correlated. Abundance was always slightly lower with the C6 than with the CANTO, except for picoeukaryotes. This could have resulted from an imperfect calibration of the analyzed volume for one of the instruments. While on the CANTO, the analyzed volume is manually calibrated, the C6 relies on a calibration every time tubing is changed following the manufacturer's recommendation. However, we have recently observed that the actual volume analyzed varies, even during the course of a day, and needs to be re-calibrated at fixed intervals using a procedure similar to the one used for the CANTO (D. Marie, unpublished data).

The slightly lower correlation coefficient observed for nanoeukaryotes could result from the difficulty to clearly distinguish the limits of the pico and nano-eukaryote populations.

In contrast, the lower sensitivity of the C6 led a drastic underestimation or even non-detection of the cyanobacteria in the upper 100 m of the water column, especially for Prochlorococcus. But the phenomenon was also present for Synechococcus although it had not been recognized previously. For both populations, corrections, such as those recommended previously [START_REF] Zubkov | Picoplanktonic community structure on an Atlantic transect from 50°N to 50°S[END_REF][START_REF] Crosbie | Abundance , distribution and flow-cytometric characterization of picophytoprokaryote populations in central (17 o S) and southern (20 o S shelf waters of the Great Barrier Reef[END_REF], did not completely solve the problem, since in some samples the corrected counts were still much lower than those obtained with the CANTO (Fig. 3C,D). The data obtained with the C6 for Prochlorococcus but also Synechococcus, should be considered with great caution even when only a part of the population is in the noise.

The oceanographic transects (Fig. 4) demonstrate that the C6 instrument would result in very serious underestimates of the cyanobacteria abundance and therefore contribution to the carbon biomass, at least in oligotrophic to mesotrophic waters. For Synechococcus, one solution could be to trigger acquisition on PE fluorescence (FL2) which is general quite strong, but this will require to run every sample twice, with the trigger on FL2 and FL3 respectively. Unfortunately this was not tested on our samples. For Prochlorococcus, the only solution to overcome the problem of low sensitivity is to increase either the excitation energy or the fluorescence detection through optical solutions [START_REF] Dusenberry | Increasing the sensitivity of a FACScan flow cytometer to study oceanic picoplankton[END_REF]Partensky et al. 1999b).

Conclusion

The increasing affordability of benchtop flow cytometers comes with limitations in some of the equipment features, such as lower detection limits. The comparison between studies of phytoplankton communities by flow cytometry should take into account the equipment used, particularly in approaches involving populations of Prochlorococcus and Synechococcus, in which low chlorophyll concentration per cell can lead to the underestimation of their abundance in the top euphotic zone. Still these benchtop flow cytometers provide reliable data for other populations such as heterotrophic bacteria and photosynthetic eukaryotes. Our study highlight the need for careful comparison between instruments before using them for large scale oceanographic surveys, using as reference the most sensitive laboratory instruments available. Table 1. Technical features of BD FACSCanto™ and BD Accuri™ C6 flow cytometers according to the manufacturer. and Synechococcus (m-x) on BD FACSCanto™ and BD Accuri™ C6. For each distribution, it is indicated whether the cells were in the noise (Noise) or whether a correction was needed (Corr.). In the depth profiles (f, l, r, x), solid symbols represent samples for which no correction was needed; grey symbols indicate samples for which we applied a correction (see Materials and Methods); samples within noise were removed. -------------------------------------------------------

Tables

Features

Fig. 1 .

 1 Fig. 1. Stations sampled in the South Atlantic Ocean off Brazil during the CARBOM cruise in 2013.Profiles: transect 1 (TR1, shaded circles) ; transect 2 (TR2, shaded triangles) and surface sampling, transect 3 (TR3, shaded squares). The grey scale on the right indicates bottom depths.

Figure 2 .

 2 Figure 2. Examples of depth profiles (St. 100 and St. 114) of normalized red fluorescence distribution (relative cell number versus chlorophyll fluorescence) and cell abundance for Prochlorococcus (a-l)

Figure 3 .

 3 Figure 3. Relationship between abundance measurements performed with BD Accuri™ C6 and BD FACSCanto™ (in cells.mL -1 ): (a) HNA bacteria; (b) LNA bacteria; (c) Prochlorococcus; (d) Synechococcus; (e) picoeukaryotes and (f) nanoeukaryotes. For Prochlorococcus and Synechococcus: 'no correction': solid circles, 'correction': grey squares; 'cells in noise': open triangles. The coefficient of determination and the equation are indicated on each graphic. The regression line calculated from 'no correction' samples is marked in black. All the slopes differed significantly from 1 (p < 0.0001), except for LNA bacteria (p = 0.049).

Figure 4 .

 4 Figure 4. Vertical abundance distribution (cells.mL -1 ) for measurements with BD FACSCanto™ (left column) and BD Accuri™ C6 (right column): Prochlorococcus (a, b, e, f) and Synechococcus (c, d, g, h). Top labels correspond to station number. Sampled points are marked as: 'no correction' (solid circles), 'correction' (grey squares), or 'cells in noise' (open triangles). Figures were drawn with the Ocean Data View software (https://odv.awi.de/).

#

  R routine to correct abundance when populations are partly in noise.R codeThe code below describes how to implement an R routine to correct the abundance of picoplanktonic populations based on their red fluorescence distribution. All libraries used here are freely available from R repositories. The input file used in this examples is named as Pro_C6.txt (See input file example File S1). This file has been created by exporting FL3 (chlorophyll) histogram from the Flowing Software (http://www.flowingsoftware.com) combining different samples into a single file. The first column contains the channel number and each following column corresponds to a different sample with rows corresponding to cell counts in each channel. Such a file could be created with any flow cytometry software. After running the cyto_plot function, a pdf output file is created named "Pro_C6.txt 1.0 .pdf" which contains all histograms from the input file (see File S3)and the file statistics (sample, uncorrected and corrected total cell abundance) are available as a data frame in the R session (see example at bottom of this file)# Example of use of cyto_plot function (run first the R code below to define the necessary functions)stats_Pro_C6<-cyto_plot("Pro_C6.txt", decades_C6, channel_min_C6, xmin_C6, xmax_C6) Set the working directory where the files are located setwd ("C:/My Documents/cytometry data/") format_format(big.mark = "", decimal.mark = ".", scientific = TRUE) # -

  determines in which case we are ("no correction", "cells in noise" or "correction")and return the corrected cell abundance in the latter case. cell_correct<-function(channel, cell_number, cell_smooth) { df<-data.frame(channel, cell_number, cell_smooth) # create a data frame i_min<-which.min(channel) # determine the minimum channel i_max<-which.max(channel) # determine the maximum channel i_cell_max<-which.max(cell_smooth) # determine in which channel is the histogram mode # "no correction" : cell abundance in the first channel is 5 times lower than abundance at the maximum of the histogram if (cell_smooth[i_cell_max]>5*cell_smooth[i_min]) {cell_correct<-"no correction"} # "cells in noise" : maximum of cell abundance is in the first channel else {if (i_cell_max==i_min) {cell_correct<-"cells in noise"} # "correction" : all the other cases, we then apply a correction by computing the total cell abundance as twice the number of cells in the channels right of the histogram maximum else {cell_correct<-2*sum(cell_number[i_cell_max:i_max]plots a set of histograms for the input samples,saves the graphics as a pdf file and compute the total cell abundance indicating whether a corrections is needed or not. It returns a dataframe containing three columns : sample, cell_tot, cell_tot_correc (see top of this file for an example) cyto_plot<-function(file_name,decades,channel_min,xmin,xmax) { channel_max = 500 # this the number of channels provided as output of the Flowing Software histo<-read.delim(file_name) histo<-histo[histo$channel>=channel_min,] histo_melt<-melt(histo, id.vars=c("channel"),variable.name = "sample", value.name = "cell_number") # smooth histogram using default R smoothing function histo_melt<-ddply(histo_melt,c("sample"), transform, cell_smooth=as.vector(smooth(cell_number))) # normalize histogram so that maximum abundance is equal to 1 histo_melt<-ddply(histo_melt,c("sample"), transform, cell_norm=cell_smooth/max(cell_smooth)) # transform log channel to linear scale for plotting if (decades==5) {histo_melt<-ddply(histo_melt,c("sample"), transform, fluo=(10^5)^(channel/channel_max))} else {histo_melt<-ddply(histo_melt,c("sample"), transform, fluo=(10^7)^(channel/channel_max))} # plots histograms using 5 columns histo_plot<-ggplot(histo_melt, aes(fluo,cell_norm)) + geom_line() + theme_bw () + facet_wrap(~ sample, nrow=21, ncol=5) + xlab("Chlorophyll")+ylab("Relative cell number") + scale_x_log10(limits=c(xmin,xmax), labels=point) # save plots as pdf ggsave(plot=histo_plot, filename=paste(file_name," 1.0 .pdf",sep=""),width = 15, height = 4, scale=2, units="cm") # compute uncorrected and corrected total cell number calling the cell_correct function defined above stats<-ddply(histo_melt,c("sample"),summarise, cell_tot=sum(cell_number),cell_tot_correc=cell_correct(channel,cell_number,cell_smooth))print(paste("# of decades:",decades,"minimum channel : ",channel_min, "xmin : ", xmin,

Table 2 .

 2 Number of samples assigned as 'correction', 'no correction' and 'cells in noise' for each picocyanobacteria group and equipment tested.

		No correction	Correction	Cells in noise
	CANTO -Prochlorococcus 74	27	1
	C6 -Prochlorococcus	27	20	55
	CANTO -Synechococcus	101	1	0
	C6 -Synechococcus	58	29	15

Acknowledgments

We thank the officers and crew of the R/V Alpha Crucis for logistical support. Financial support for this work was provided by the CNRS Groupement de Recherche International (GDRI) "Diversity, Evolution and Biotechnology of Marine Algae", the COFECUB-CAPES project "Pico Brás" (TE 871-15) and FAPESP -Fundação de Amparo à Pesquisa do Estado de São Paulo (2012/04800-9 and 2014/15242-2).

Supplementary material is available at https://figshare.com/s/a9499d9ab4f4740eb576 .

Supplementary figure legend