Romain Noël
email: romain.noel@univ-eiffel.fr

Laurent Navarro
email: navarro@emse.fr

Guy Courbebaisse
email: guy.courbebaisse@insa-lyon.fr

Lattice Boltzmann Method for Mathematical Morphology: Application to Porous Media

Keywords: Lattice Boltzmann Method, Mathematical Morphology, Image Processing, Porous Media, Porosity, Permeability

The LBM (Lattice Boltzmann Method) is often used in CFD (Computational Fluid Dynamics) for efficient fluid flow simulations. Computation of the permeability of a porous media from direct simulations is a common application which benefits from the ability of the LBM (Lattice Boltzmann Method) to embed porosity parameters. The MM (Mathematical Morphology) is widely used in image processing as the theoretical aspects guaranty robust algorithms for geometrical characterization of shapes appearing in images. The MM is commonly used to compute porosity from porous media images. The union of these two methods has been recently done through the LB3M (Lattice Boltzmann Method for Mathematical Morphology).

The present work extends the LB3M to the extraction of porosity and pores segmentation from images. In order to benefit from the full capacity of the LB3M, it is necessary to reformulate and adjust the algorithms in a new paradigm. Thus, the underlying concept and algorithms required for computing the different previous information are detailed. Moreover, a comparison is provided between the permeability resulting from the CFD and MM both implemented by using the LBM.

To sum up, this work emphasizes the full capacity of the LB3M to obtain complex transformations and operations issued from the MM theory through completely new and innovative algorithms. The herein challenge is to highlight the abilities of the LB3M to match with physical phenomenons. Indeed, the LB3M keeps the advantages from the MM such as a complete theory, fast convergence, scalability, robustness, etc. while adding the power of the LBM: statistical physics origins, partial differential equation solver, intrinsic properties of parallelization, efficiency, etc.

INTRODUCTION

The LBM * is now widely used for computational fluid dynamics, because of its advantages in terms of parallel computing facilitation, competitive computation times, natural micro to macro approach, physically interpretable simulation,... Since some decades, image processing benefits several of the advantages of the LBM as a numerical method tool. However, until recent works, the connection of image processing and physics in the LBM framework did not exist.

The LBM is a method that uses particle density on a regular network to solve fluid dynamics equations. [START_REF] Succi | The Lattice Boltzmann Equation: For Complex States of Flowing Matter[END_REF] The network is fixed and the description of phenomena is Eulerian by nature. Thus, the LBM network can be seen as an image in D dimensions where grey-levels correspond to density. This opens the way to the definition of many operators and image processing methods in the LBM framework.

Traditionally, image processing using the LBM was more oriented to methods that need somehow a diffusion process. Indeed, the LBM with BGK † approximation can produce a linearized diffusion process. Thus, various algorithms based on diffusion and anisotropic diffusion have been ported in the LBM framework. [START_REF] Jawerth | Lattice Boltzmann models for anisotropic diffusion of images[END_REF][START_REF] Chen | An Anisotropic Diffusion Model for Medical Image Smoothing by Using the Lattice Boltzmann Method[END_REF] Applications like denoising, inpainting, or segmentation have been successfully developed. However, one of the most important tools of image processing is the MM. ‡ MM is used in nearly any image processing algorithm, even for a small pre-or post-treatment. Most of image processing toolboxes of any computing language embed MM tools.

The LB3M § has been proposed in 2019, [START_REF] Noël | Lattice Boltzmann Method & Mathematical Morphology[END_REF][START_REF] Noël | The Lattice Boltzmann Method for Numerical Simulation of Continuum Medium Aiming Image-Based Diagnostics[END_REF] in a paper defining the basic operators and demonstrating the feasibility of using LB3M for bone tissue remodeling simulations. The principle of the proposed work was to define basic MM operators using a modification of the equilibrium distribution in the BGK approximation. This led to dilation and erosion operators defined using pixel values as particle density.

The main advantage of using LB3M is to do any operation in a LBM framework when image processing and numerical simulations are required. An example of a complete pipeline using the LBM can be found in other works. [START_REF] Noël | Lattice Boltzmann method for modelling of biological phenomena[END_REF] Moreover, this gives the opportunity to define in-between operations and specific operators that can lead to "physical" image processing. In the example of bone tissue remodeling, the idea was to use the dilation operator in a controlled manner to increase bone volume in scanner images. By controlling the dilation operator, it is possible to embed some particular rules of bone biology in the process.

Porosity is defined simply as the proportion of pores, or the absence of considered matter, in a solid material. It is a physical quantity computed as the ratio between volume of pores (holes) and total volume of a material, comprised between 0 and 1. The computation of porosity is important, as the pores, if not closed, can be filled with liquid, smaller particles, biological cells,... Thus, the computation of mass density of a material, for example, is highly impacted by porosity.

In this paper, we propose to compute the porosity of materials using LB3M on material images. We first show how the porosity can be calculated using mathematical morphology operators based on the morphological distance transformation on a synthetic example. Then, we present a method to compute fluid permeability using the LBM, and we detail the different algorithms of the LB3M operators. We also propose a validation of the process via the computation of porosity and permeability on the synthetic example and a real case of clay spheres from a heat exchanger. This new approach for the computation of the porosity allows the mixing of numerical modeling with image processing, and it gives the opportunity to compute porosity and permeability in one single framework.

The paper is organized as follows. Section 2 presents the basics of mathematical morphology and two methods for computing porosity (SKIZ ¶ and watershed), and the basics of LBM for fluid permeability computation. The implementation in the LB3M framework is also detailed. Section 3 is dedicated to results on a synthetic case and the application to porous medium images. Results are discussed in section 4, and conclusions and perspectives are given in section 5.

METHOD

Mathematical Morphology Porosity

There exists several methods to compute the porosity based on images of a material. 7 Among all of these methods, those based on mathematical morphology [START_REF] Matheron | Elements pour une theorie des milieux poreux[END_REF][START_REF] Najman | Mathematical Morphology: From Theory to Applications[END_REF] rely on pure geometrical characteristics, which is very important for the profound understanding of the porosity. This section first recalls the basics of the MM, and † Bhatnagar, Gross and Krook ‡ Mathematical Morphology § Lattice Boltzmann Method for Mathematical Morphology ¶ SKeletonization by Influence Zones then details two reliable methods, SKIZ transformation and watershed using markers, to perform a segmentation allowing the computation of the porosity ϕ with the ratio:

ϕ = V pores V total . (1)
A grey-level image ρ in D dimensions, is mathematically defined as a function of R D in R. Furthermore, a greyscale structuring element g is defined as a function that goes from R D into R and such as the

set x ∈ R D | g(x) = {∞} ∪ {-∞} is bounded.
Then the morphological dilation of a grey-level image ρ by a functional structuring element g, is denoted (ρ⊕ g) and is defined by the relation

(ρ⊕ g) (x) = sup y (ρ (y) + g (x -y)) . (2)
Oppositely, a morphological erosion of a grey-level image ρ by a functional structuring element g, is denoted (ρ g) and is defined such as

(ρ g) (x) = inf y (ρ (y) -g (y -x)) . (3)
Examples of these two operations applied on a synthetic binary image (see fig. 1) can be seen in fig. 2. It is worth mentioning that all figs. 2 to 7 used to illustrate the MM tools were obtained by using LB3M operations only. From these two basic operators, more complex operators and a whole mathematical theory can be derived. More specifically, it is possible to define a bounded dilation for the reconstruction of the function h, using the reference image ρ and the functional structuring element g:

R (h, g, ρ) = [min {(h⊕ g) , ρ}] ∞ . (4)
One can also, using the erosion operator on a binary image ρ, compute a map representing the geodesic distances for all the shapes in the image:

d (ρ) (x) = ( ρ b) ∞ (x) , ( 5
)
where b is the functional structuring element associated to the geodesic ball (i.e. related to the chosen norm) and ρ is the rescaled binary image going from 0 to +∞. This operation can be seen as an iterative infinite erosion. The distance transformation image from the synthetic image of fig. 1 is shown fig. 3. False colors are used in this image, blue represents low values and yellow represents high values. Thus, the more yellowish the pixels are, the thicker the image's shapes are. Starting from the distance transformation, a watershed operator 10 can be defined. It consists in the progressive flooding of the distance map with virtual water, until the different regional minima of the map overflow. The watershed gives a segmentation of the image based on barriers constructed when the merging of two or more water flooding sources occurs.

We first define the so-called Influence Zones (IZ), which are the sets of points that are closer to each of the tags T :

IZ n (ρ, T n )(x) = n if d (T n , ρ) (x) < d (T k , ρ) (x) ∀ k 0 else . (6)
Then, the watershed of the distance map is computed as: One solution to cope with the over-segmentation drawback is to perform the watershed segmentation from markers. Several markers types can be used, but a straightforward choice is the maxima of the height of the distance map given by:

lim n→∞ W n (ρ) = lim n→∞ IZ n (ρ = n, W n-1 ) + ([ρ = n] -R (ρ = n -1, b, ρ = n)). (7
Maxima (ρ, h) (x) = ρ(x) -R (ρ -h, b, ρ) (x). ( 8
)
This allows to drive the watershed using the information of the image itself. The markers obtained from the distance map (fig. 3 Then, integrating the markers to the watershed operator gives:

lim n→∞ W n (ρ) = lim n→∞ max{IZ n (ρ = n, W n-1 ) , T }. ( 9
)
Figure 6 shows the results of the segmentation using watershed with markers. The use of the markers clearly reduces the over-segmentation drawback. It even solves the problem entirely in our example. Another method that can be used to do the segmentation is the SKIZ. It is defined as the boundary of all IZ starting from the tags T :

SKIZ(ρ, T ) = ρ - n 1 n IZ n . (10)
An example of the segmentation obtained using SKIZ transformation is shown fig. 7. One can see differences between watershed with markers and SKIZ transformation. Indeed, in the SKIZ transformation, the lines connecting the circle are not shortest paths between them, i.e. these lines are not aligned to the centers of the circles. This is due to the fact that SKIZ transformations found separation lines in the middle of geodesic distance between the tags without considering the topology brought by the distance map.

LBM: Fluid Permeability

The LBM is a method for numerical analysis of physical processes and phenomena. The LBM is an alternative to classical differentiation schemes allowing a direct numerical implementation of a virtual microscopic physics and an a posteriori averaging for estimating macroscopic parameters. Considering the Boltzmann equation,

∂ f ∂t + ξ • ∇ x (f ) = Ω (f , f ) , (11)
the collision operator Ω (f , f ) is expressed as a complex non-linear integral. To compensate this disadvantage a common way consists in linearizing the collision term around its local equilibrium solution. Ω (f , f ) is often replaced by the so-called BGK collision operator. [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] After introducing the BGK approximation, the Boltzmann equation can be written in a completely discretized equation providing the lattice Boltzmann equation in a BGK formulation, called the LBGKE:

‖ f i (x + ξ i ∆t, t + ∆t) = f i (x, t) -∆tω f i (x, t) -f (eq) i (x, t) , (12)
where ∆t is a time increment and ξ i ∆t a space step, f i is the particle velocity distribution along the i-th direction and ω is the inverse of the relaxation time (or relaxation frequency). In addition, f (eq) i

is the equilibrium distribution function computed as:

f (eq) i (ρ, v) = w i ρ 1 + ξ i • v c s 2 + (ξ i • v) 2 2c s 4 - (v) 2 2c s 2 , ( 13
)
in which v is a macroscopic velocity of particles in a medium, ξ i are velocity vectors in the specific lattice model, ρ is the macroscopic mass density, and c s represents speed of sound. Another specificity of the method is the fact that the LBM is usually solved in two sub-steps: a collision step and a streaming step expressed successively as follows:

f i c (x, t) = f i s (x, t) -∆tω f i s (x, t) -f (eq) i (x, t) , (14)
f i s (x + ξ i ∆t, t + ∆t) = f i c (x, t). (15)
‖ Lattice Boltzmann-BGK Equation

The connections of the distribution functions to macroscopic quantities for the density, the momentum flux and the stress tensor are defined as follows:

q i=0 f i = ρ, (16)
q i=0 ξ i f i = ρv, (17)
q i=0 (ξ i -v) ⊗ (ξ i -v)f i = Π = p (0) -1 + ∆tω 2 σ (1) , (18)
where p (0) and σ [START_REF] Succi | The Lattice Boltzmann Equation: For Complex States of Flowing Matter[END_REF] are given by the Chapman-Enskog [START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction, and Diffusion in Gases[END_REF] and are equal to the ideal pressure and ∇ x (v), which allows to relate 1 + ∆tω 2 to kinematic viscosity. Also note that all the quantities given in this paper are using the LBM formalism and therefore are dimensionless. As we want to estimate the permeability of material, the passage of microscopic flows at the pore scale is governed by NSE:

* * ∂ t (ρv) + ∇ x • (ρv ⊗ v + p -ν∇ x v) = 0. ( 19
)
In the same time, it is well established that LBM is an adapted technique for numerical simulations particularly in the field of fluid mechanics. In this framework, the LBM in BGK approximation is implemented for solving the NSE as the NSE can be recovered from a Chapman-Enskog expansion. [START_REF] Malaspinas | Simulation of a generalized Newtonian fluid by the Lattice Boltzmann Method[END_REF] In addition in the context of porosity estimate, the NSE has to be adapted for capturing the fluid flow inside a porous medium. Indeed, the fluid in a porous medium is not flowing freely due to its interaction with the solid matrix, which differs from the continuum hypothesis of the NSE. A common adaptation is obtained through Darcy's law, which reads:

q = - K ρν (∇ (p) -ρg), (20)
where q is Darcy's speed i.e. the volume flow of the fluid, p is the pressure, ρ is the mass density, ν is the kinematic viscosity of the fluid, g is the gravity acceleration field and K is the permeability characterizing the porous medium.

The permeability of a porous medium and its porosity are intrinsically correlated, but are not identical. The Kozeny-Carman equation, [START_REF] Hommel | Porosity-Permeability Relations for Evolving Pore Space: A Review with a Focus on (Bio-)geochemically Altered Porous Media[END_REF] is a simple expression of the relationship between the porosity and permeability. This latter can be expressed as:

K = Φ 2 s D 2 p 180 ϕ 3 (1 -ϕ) 2 , ( 21
)
where is Φ s the sphericity and D p is the characteristic diameter of the particles.

LB3M Algorithmic

This section presents the basics of the LB3M. Starting from the classical mathematical morphology operator definitions, we show that these operators can be derived from the LBM equation. To do this, we consider an equivalence between the image and the density in the LBM framework, both denoted ρ.

On the one hand, the MM is partly based on dilation and erosion operators and relies on lattice theory. On the other hand, the LBM is by definition based on lattice theory and a recent work [START_REF] Noël | Lattice Boltzmann Method & Mathematical Morphology[END_REF] has proved the ability of the LBM to reproduce exactly the morphological dilation and erosion, leading to the LB3M. Therefore, the LB3M is able to reproduce all the operators of the MM within the LBM formalism. This formalism brings an adaptation of the algorithms, and these adaptations to the morphological tools mentioned in the previous section are detailed in this section.

As previous work mentioned, [START_REF] Noël | Lattice Boltzmann Method & Mathematical Morphology[END_REF] the ability of the LB3M to reproduce dilation relies on the modifications of the equilibrium distribution and the norm used to compute the macroscopic moments. In other words, by replacing the equations ( 13) and ( 16) respectively, by the following equations:

ρ (x, t) = max i {f i (x, t)} (22a) f (eq) i (x, t) = ρ (x, t) + g (ξ i ) , (22b)
the LB3M reproduces exactly (from a mathematical point of view) the morphological greyscale dilation. These modifications lead to the algorithm 1. It is worth noting the fact that parameter ω is set to 1 in order to obtain

f i , ρ dilate = ρ, ω = 1; for n ← 1 to n max do for i ← 0 to q do f (eq) i (x) = ρ dilate (x) + g(ξ i ) ; // Compute equilibrium distribution f i (x) = f i (x) -ω(f i (x) -f (eq) i (x)) ; // Collide f i (x) = f i (x -ξ i ) ; // Stream end ρ dilate (x) = max i {f i (x)} ;
// Compute the macroscopic dilation end a complete dilation, but a value in [0, 1] will lead to a partial dilation and the 0 value is resulting no dilation at all. This property is a natural generalization of the dilation operator brought by the LBM formalism and will be used later to take full advantage of the LB3M.

The erosion is the complementary operation (but not inverse) of the dilation. Very similarly to dilation, it has been proved that by replacing eqs. ( 13) and ( 16) by the equations:

ρ (x, t) = min i {f i (x, t)} , (23a)
f (eq) i (x, t) =ρ (x, t) -g (-ξ i ) , ( 23b
)
the LB3M is reproducing the greyscale morphological erosion. The pseudo-code reflecting this operator is given in algorithm 2. The same remark about the generality of the parameter ω can be made.

Algorithm 2: Erosion in LB3M

Input: image ρ, structuring element g, number of erosion n max Result: Morphological erosion Initialize variables:

f i , ρ erosion = ρ, ω = 1; for n ← 1 to n max do for i ← 0 to q do f (eq) i (x) = ρ erosion (x) -g(-ξ i ) ; // Compute equilibrium distribution f i (x) = f i (x) -ω(f i (x) -f (eq) i (x)) ; // Collide f i (x) = f i (x -ξ i ) ; // Stream end ρ erosion (x) = min i {f i (x)} ;
// Compute the macroscopic erosion end

The reconstruction transformation is mathematically close to the dilation. Thus, the algorithm 3 is close to the LB3M dilation algorithm, with some details differing. The converging character of the reconstruction is brought by a while loop, and the reference set is introduced for comparison with a min operator in the macroscopic computation.

Algorithm 3: Reconstruction in LB3M

Input: image h, structuring element g, reference image ρ Result: Morphological reconstruction Initialize variables:

f i , h 0 reconstructed = h, ω = 1; while h t reconstructed = h t-1 reconstructed do t = t + 1 ; // Increment for i ← 0 to q do f (eq) i (x) = h t-1 reconstructed (x) + g(ξ i ) ; // Compute equilibrium distribution f i (x) = f i (x) -ω(f i (x) -f i (eq) (x)) ; // Collide f i (x) = f i (x -ξ i ) ; // Stream end h t reconstructed (x) = min max i {f i (x)} , ρ(x) ;
// Compute the macroscopic reconstruction end

The algorithm 4 giving the distance map algorithmic in the LB3M formalism is even closer to the LB3M erosion, since it is an erosion with a geodesic ball as structuring element until convergence. So, without surprise, the algorithm 4 is differing from the LB3M erosion only from the while loop checking the convergence of the distance map. // Increment for i ← 0 to q do f (eq) i

(x) = d t-1 (x) -b(-ξ i ) ; // Compute equilibrium distribution f i (x) = f i (x) -ω(f i (x) -f (eq) i (x)) ; // Collide f i (x) = f i (x -ξ i ) ; // Stream end d t (x) = min i {f i (x)} ;
// Compute the macroscopic distance end

Starting from the distance map, maxima of such map can be extracted easily using reconstruction as eq. ( 8) emphasize. Such maxima extraction can be done with a simple MM formalism, even if a specific development in the LB3M could be interesting, it does not bring new concepts. This is why, no specific LB3M has been developed for this purpose in this paper. However, these maxima are essential in the present work since they are capturing tags. Then, these tags are themselves used as starting points of the IZ and watershed algorithms.

The geodesic IZ is a key algorithm for both SKIZ and watershed transformations. Interestingly, it is also the less straightforward application of the LB3M formalism. It uses both dilation and erosion operators in a single iteration. The idea is to construct step by step the geodesic distance map of all the tags in the image indistinctly of the tag value. And to discriminate to which tag belongs the minimal distance value, a conditional reconstruction of tags in the image is linked to direction from which comes the minimum. This conditional reconstruction is made using the generalization brought by the parameter ω. Thus, algorithm 5 is summarizing the strategy adopted to compute the geodesic IZ in the LB3M framework. Since a distance map is used, the algorithm requires a first distribution f i with a constant parameter ω equal to 1 and a geodesic ball b as structuring element. But since the tags are reconstructed simultaneously, the algorithm also requires a second distribution g i , a flat structuring element noted b (equal to 0 inside the geodesic ball and -∞ outside). This second distribution needs a parameter ω i depending of space and directions i to control the dilation of tags inside the reconstruction.

Algorithm 5: IZ in LB3M

Input: image rescaled ρ, tags T Result: geodesic IZ transformation Initialize variables: f i , geodesic ball b, ω = 1, distance map d 0 = ρ, g i , flat structuring ball b, ω i = 0,

IZ 0 = T ; while d t = d t-1 do t = t + 1 ; // Increment for i ← 0 to q do f (eq) i (x) = d t-1 (x) -b(-ξ i ) ;
// Compute equilibrium distribution for erosion The geodesic IZ and by extension the SKIZ are very powerful tools, but in the case of pore segmentation they are missing the topology given by the solid matrix. Or more precisely, the topology missing in the IZ is given by the distance map related to the solid domain of the porous medium. To overcome this issue, the watershed is combining the IZ with the distance map information. However, the simple watershed by flooding is using automatically all the minima from the distance map to construct the catchment basins. This automatic use of minima leads to an unwanted over-segmentation. Therefore, to workaround this disadvantage, the usual watershed by markers is an appropriate tool. As for the maxima, the watershed algorithm adopted in this paper is more an application of the MM equation than a reformulation in the LB3M paradigm. But unlikely to the maxima, the watershed algorithm deserves a bit more attention since it is the final results we obtained here, and is presented in algorithm 6. for n ← 0 to max(ρ) do W n (x) = max IZ(ρ = n, W n-1 )(x) , T (x) end

f i (x) = f i (x) -ω(f i (x) -f (eq) i (x)) ; // Collide f i (x) = f i (x -ξ i ) ; // Stream end d t (x) = min i {f i (x)} ; // Compute the macroscopic distance for i ← 0 to q do g (eq) i (x) = IZ t-1 (x) + b(+ξ i ) ; // Compute equilibrium distribution for dilation ω i (x + ξ i ) = 1 if f i (x) = d t (x) 0 otherwise ; // Compute ωi according to the contribution of f i to d t g i (x) = g i (x) -ω i (x)(g i (x) -g (eq) i (x)) ; // Collide g i (x) = g i (x -ξ i ) ;

RESULTS

Comparison between LBM and LB3M

A conversion is needed to compare the geometrical results from the LB3M and fluid results from the LBM. This common quantity that let us compare geometrical and fluid results is the permeability. Indeed, from the porosity computed from the LB3M, permeability can be computed thanks to the Kozeny-Carman equation, see eq. ( 21) ; while the permeability can be estimated from a fluid flow through a porous medium thanks to Darcy's law, see eq. ( 20). In the sake of consistency and simplicity, the previous example image of a porous medium (see fig. 1) is used for computations with both LBM and LB3M.

Starting from the LB3M watershed segmentation obtained in fig. 6, the area of each pore is extracted by counting the number of pixels having the tag label. Then, all the pore areas are summed to obtain the total pore area that is compared to the total image, to compute the porosity thought the eq. ( 1). Thus, we computed a total pore area of 20 715 pixels for total image area of 33 000 pixels. This leads to a porosity of 0.627 72. Therefore, the Kozeny-Carman equation yields the permeability of 10.99. On the LBM side, the image is extruded in the third dimension to obtain the 3D geometry of the porous medium. Thus, the beads and walls on the image are set as rigid impermeable boundaries. A very common in LBM formulation of boundary condition is used to model these walls: a full-way bounce-back condition. The bottom and top images are used to impose extra pressure boundary conditions. Therefore, in the pressure gradient is controlled by the boundary conditions, such as the viscosity. Then, we can measure the volume flow of the fluid by integrating the speed magnitude (see fig. 8) over a cross section. Therefore, using the Darcy's law, we can estimate the permeability of the previous porous medium image around K = 10.33.

The permeability values computed from the LB3M and LBM are close without being exactly equal. This proximity between the values indicates a good fidelity of the results obtained by both LB3M and LBM. Indeed, the difference between these values can have multiple origins like the estimations of the parameters, different sensibility to the parameters, incertitude etc. These values are, in our case, more indicators than ground truth.

Application to Clay Spheres

In this section we propose the application of the porosity calculation to clay material intended to be integrated to the design of a heat exchanger for data center cooling. The system is a direct contact heat exchanger built as a vertical packed bed, composed of earthen spheres. Details can be found in Bouchenna's paper. [START_REF] Bouchenna | Heat exchanger design based on earthen materials[END_REF] Images used here correspond to one layer of spheres of three different sizes as shown fig. 9. The spheres are disposed in a cylindrical container. Computation of porosity is of main importance in this application, in order to characterize the energy buffering capacity of the heat exchanger. The methods described in section 2 are based on binary images. Thus, the first step is to segment the color images of the spheres. To do so, we transformed the images from RGB color space to L * a * b * color space and we kept the b * channel only. Indeed, the blue-yellow opposition is well adapted to the discrimination of the color of the spheres. Then, we performed a classical Otsu automatic thresholding [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF] From these binary images, it is possible to calculate the distance map with eq. ( 5). The resulting distance maps from the three images of corresponding to the three sizes of spheres are shown in fig. 11. One can observe that there are differences in amplitudes. These differences are due to the arrangements of the spheres: one can see that if the arrangement fits well in the container, the spheres present symmetrical patterns like in fig. 10c. A summary of all the results obtained with the watershed by markers (see algorithm 6) for the three different sizes of spheres is proposed in table 1. The cylinder area (diameter of the container), the pore areas sum, the number of pores, the min-max of pore areas, the mean value and standard deviation, and the global porosity are computed.

All the values in the table 1 except cylinder area and porosity have to be considered regarding the cylinder area values, as the images do not have exactly the same resolution due to the operator taking the picture. One can observe that global porosity is increasing with spheres diameters, which is expected. However, the other values (pore areas sum, number of pores, Min-Max of pores areas, mean value, and standard deviation) allow a finer description of the material's geometric properties linked to porosity than the global porosity alone.

DISCUSSIONS

This paper proposes to apply the LBM to the computation of porous media characteristics: porosity with the LB3M and permeability with LBM fluid simulations. The mathematical morphology is widely used to compute porosity from images of materials. More specifically, the SKIZ transformation and the watershed using markers are good candidates, as they both avoid the over-segmentation phenomenon. We recalled these methods while illustrating the different steps on a synthetic example image. We also detailed how to compute fluid permeability using the LBM. Then, we proposed several LB3M algorithms to perform all the different steps of SKIZ transform and watershed using markers.

To illustrate the method we proposed to apply the porosity calculations computations on clay spheres that make parts of a heat exchanger designated for the cooling of data centers. This application has been chosen since the determination of porosity is a major parameter allowing the good design of such heat exchanger. In a previous work 15 a limitation appears from the incertitude in the measurement method used to estimate the porosity. Here, we offer to tackle this limitation with an original manner.

The main advantage of the LB3M is the use of a single framework for all the porosity characteristic computation. A previous work [START_REF] Noël | Lattice Boltzmann Method & Mathematical Morphology[END_REF] shows the interest of a single framework with an application on bone tissue growth from scanner images. Indeed, even if all the steps could be realized using classical mathematical morphology, it gives the opportunity to mix fluid simulation with image analysis. This is a really important characteristic of the method, as if the spheres are moving during fluid simulations for example, one can observe the variation of porosity without the need to export images to analyze them in another framework. This is particularly true for heat exchangers which require thermal flow simulations. A kind of simulation for which the LBM is a method of choice.

Of course, the heat exchanger evoked previously is a 3D stack of clay beads. So, a more realistic application and computation should be in 3D. However, the acquisition of such image of the inner stack requires more advance devices (CT-scan, MRI, etc.). But the comparison with 3D LBM fluid flow would also be more detailed with such 3D images. We hope to have access to this kind of data in a close future to pursuit further comparisons between LBM and LB3M.

Another application of the watershed that worth being mentioned, since related to fluid flow, is the extraction of porosity and geometries to construct a pore network. Such networks are very useful to build a PNM. † † The PNM is a numerical method alternative to LBM, but is rather limited to simple fluid flow, since its accuracy is limited by some important simplifications. On the other hand, these simplifications make the PNM a very fast method. The PNM requires to know the volume, connectivity and geometry of the pores but also of the throats: information often obtained thanks to MM. Sometimes, the PNM is even sometimes calibrated by using a LBM simulation, showing again the possibility of the LB3M at the frontier between LBM and MM.

Concerning the computation of SKIZ, the LB3M really shows its powerfulness by using the conditional dilation as shown in algorithm 5, they take full advantage of the LB3M formalism and generality. Indeed, a straightforward approach could have been to compute a different geodesic distance map for each tag, store them, then affect to the point IZ(x) the tag label associated to map with the minimum distance value at x. However, this approach would have been more computationally expensive since the distance maps would have been computed very far from the tag zone (unnecessary area if other tags are closer); plus, it would have been more memory demanding since storing as many distance map as the tag number is much larger than only one common distance map. Thus, such conditional dilation approach is more natural with the LB3M. Even if, it is possible with the traditional MM, it would remain less adapted.

The watershed computation also benefits from the previous advantage, since it uses the IZ computations. Nevertheless, a remark can be made: the algorithm 6 presented and used here is not optimal. Indeed, since the watershed algorithm uses many times the algorithm IZ, it initializes the variable of IZ very time, while they are close to the ones used in the previous iteration. This is far from being optimal because the allocation/deallocation of these variables is time consuming with any purpose. Therefore, a first improvement would be to incorporate the IZ procedure directly into the watershed algorithm that reuses the variable from their previous state rather than reallocate and reinitialize them. A second improvement could be to include topological levels directly into the geodesic distance rather than looping over levels, but this requires more precautions.

CONCLUSIONS

The LBM is a well-adapted method to solve Navier-Stokes Equations. However in the last decade LBM has been redefined as a potential method for the image processing domain. For instance, denoising and segmentation implementation have been performed to treat images sequences. As the LBM is intrinsically highly parallelizable, the calculation times are significantly reduced. This paper reminds recent works when applying LBM for implementing mathematical morphology operators. In the continuity, this paper presents new investigations performed such as watershed segmentation by markers and SKIZ segmentation applied to the estimation of the porosity of materials. To conclude, the originality of this paper lies in the comparison of the results provided through the solving of the NSE and through the image processing operators, the two approaches having been implemented through the LBM. The results are fruitful and of course, need to be confirmed on different applications notably concerning the gain on the computation time.

Figure 1 :

 1 Figure 1: Original image of an example porous medium.

  (a) Dilation using a flat square structuring element (b) Erosion using a Gaussian structuring element

Figure 2 :

 2 Figure 2: Dilation and erosion of the example porous medium.

Figure 3 :

 3 Distance transformation of the porous medium.

)

  One might note that the part ([ρ = n] -R (ρ = n -1, b, ρ = n)) is the new minima at level n of the function ρ. An example of the watershed operation performed on the original image of fig. 1 is shown fig. 4. This operation is automatic, thus sensitive to noise and digitization parameters. One can observe the over-segmentation drawback, characterized by too numerous small regions.

Figure 4 :

 4 Watershed segmentation by flooding of the porous medium, leading to an over-segmentation.

Figure 5 :

 5 Markers obtained from maxima of the distance transformation.

Figure 6 :

 6 Watershed segmentation by flooding of the porous medium from markers.

Figure 7 :

 7 Segmentation by SKIZ transformation of the porous medium.

Algorithm 1 :

 1 Dilation in LB3M Input: image ρ, structuring element g, number of dilation n max Result: Morphological dilation Initialize variables:

Algorithm 4 :

 4 Distance map in LB3M Input: image rescaled ρ Result: Morphological distance map Initialize variables: f i , geodesic ball b, ω = 1, distance map d 0 = ρ; while d t = d t-1 do t = t + 1 ;

  end IZ t (x) = min max i {g i (x)} , ρ ; // Compute the macroscopic reconstruction of tags end

Algorithm 6 :

 6 Watershed by markers in MM Input: image ρ, tags T Result: Watershed transformation Initialize variables: W -1 = T , ;

Figure 8 :

 8 Figure 8: Color map slices of the speed magnitude from 3D LBM fluid flow through a porous media.

Figure 9 :

 9 Figure 9: Original color RGB images of a porous medium.

  and used a mask to keep only the relevant part of the image, i.e. inside the container. The resulting binary images are shown fig. 10.

Figure 10 :

 10 Figure 10: Binarisation of the original images with an Otsu thresholding.

Figure 11 :

 11 Figure 11: Distance transformations of the binary images.

Figure 12 :

 12 Figure 12: Segmentation by SKIZ transformations of the distance maps, starting from the markers.

Figure 13 :

 13 Figure 13: Watershed segmentation by markers of the distance maps.

Table 1 :

 1 Max of pores areas 1803 -43 194 15 709 -144 637 19 550 -183 189 Porosity computed from the watershed.

		D15	D30	D45
	Cylinder area	4 582 540	4 677 345	3 382 183
	Pore areas sum	1 005 783	1 337 418	984 995
	Number of pores	107	27	12
	Min-Mean value	9399.84	53 496.72	82 082.92
	Standard deviation	7578.37	28 007.32	59 282.62
	Global porosity	0.2195	0.2859	0.2912

* * Navier-Stokes Equations

† † Pore Network Model

ACKNOWLEDGMENTS

This work has been partially funded by the French National Research Agency (Agence Nationale de la Recherche) via the LBSMI project ANR-15-CE19-0002. This support is gratefully acknowledged.


P G Albert

B Giaccio

R Isaia

A Costa

E M Niespolo

S Nomade

A Pereira

P R Renne

A Hinchliffe

D F Mark

PG Albert

V C Smith

Affiliations

Evidence for a large-magnitude eruption from Campi Flegrei caldera

published or not. The documents may come

L'archive ouverte

INTRODUCTION

Near-vent volcanic successions often provide fragmentary records of past explosive volcanism due to the intense volcano-tectonic and syn-eruptive sedimentary processes they experience. This means they are not always representative of the full range in scale and tempo of past activity. Statistical analysis of existing eruptive databases verifies underrecording, including large magnitude eruptions [START_REF] Kiyosugi | How many explosive eruptions are missing from the geological record? Analysis of the Quaternary record of large magnitude explosive eruptions in Japan[END_REF], which is problematic for hazard assessments. Fortunately, explosive eruption deposits are recorded as ash (tephra) layers in sedimentary records. In the Mediterranean region the distal tephrostratigraphic record has been widely used as a key tool for studying past volcanism, reconstructing long-term eruption histories and constraining the tempo and dynamics of explosive eruptions (e.g., [START_REF] Paterne | Explosive activity of the South Italian volcanoes during the past 80,000 years as determined by marine tephrochronology[END_REF][START_REF] Wulf | Tephrochronology of the 100 ka lacustrine sediment record of Lago Grande di Monticchio (southern Italy)[END_REF].

Although some of the most widespread Mediterranean ash layers are correlated to nearsource deposits, e.g., the Y-5 marine tephra linked to the Campanian Ignimbrite (CI) eruption of Campi Flegrei caldera (CFc) at 40 ka (e.g., [START_REF] Giaccio | High-precision 14 C dating and 40 Ar/ 39 Ar dating of the Campanian Ignimbrite (Y-5) reconciles the time-scales of climaticcultural processes at 40 ka[END_REF], many other widely dispersed tephra do not have counterparts near source. One of the most widespread layers not confidently correlated to a particular eruption deposit is the Y-3 tephra. This lastglacial tephra layer was first reported in the sediments of the Ionian Sea [START_REF] Keller | Explosive volcanic activity in the Mediterranean over the past 200,000 yrs as recorded in deep-sea sediments[END_REF], and has been traced across the Mediterranean (Fig. 1A; Table S1). A precise age of 28,690-29,420 cal yrs BP was obtained for the Y-3 in the Tenaghi Philippon peat sequence, Greece [START_REF] Albert | The Y-3 tephrostratigraphic marker revisited: new diagnostic glass geochemistry, improved chronology and climatostratigraphic interpretations[END_REF]. The glass compositions of the Y-3 tephra are consistent with CFc, but within the so-called 'Tufi Biancastri' (TB) succession [START_REF] Orsi | The restless, resurgent, Campi Flegrei nested caldera (Italy): Constraints on its evolution and configuration[END_REF][START_REF] Tomlinson | Geochemistry of the Phlegraean Fields (Italy) proximal sources for major Mediterranean tephras: implications for the dispersal of Plinian and co-ignimbritic components of explosive eruptions: Geochim[END_REF] encompassing all the eruption deposits between the ~15 ka Neapolitan Yellow Tuff (NYT; [START_REF] Deino | The age of the Neapolitan Yellow Tuff caldera-forming eruption (Campi Flegrei caldera -Italy) assessed by 40 Ar/ 39 Ar dating method[END_REF] and CI (Fig. 1B). No deposit with the full range of compositions (phono-trachytic to trachytic) diagnostic of the distal layer had been found at source (see [START_REF] Albert | The Y-3 tephrostratigraphic marker revisited: new diagnostic glass geochemistry, improved chronology and climatostratigraphic interpretations[END_REF]. Here we report the characteristics of a mid-proximal TB deposit that correlates with the Y-3 distal layer, and helps verify the eruption was a large magnitude event.

PROXIMAL STRATIGRAPHY AND CHEMICAL COMPOSITIONS OF THE DEPOSITS

The distal Y-3 tephra has two volcanic glass compositional end-members, higher-SiO 2 trachytic (Component 1) and lower-SiO 2 phono-trachytic (Component 2) glasses, with the phono-trachytic glasses being particularly distinctive. Chemical analyses (methods see Data Repository 1) of TB units identified in a borehole (S-19) and an outcrop at Ponti Rossi (PR; Fig. S1), northeast of CFc in central Naples, show the two compositional end-members are consistent with the distal tephra (Fig. 2; Table S2).

The S-19 TB sequence is 14 meters thick with 6 primary tephra units that range from a few cm to 3.4 m in thickness. The thickest unit is the third above the CI, found at 25. 90-22.50 m (3.40 m thick) and it presents the glass compositions of the distal Y-3 tephra (Fig. 2). The TB sequence at PR has eight eruption units separated by eruption hiatuses marked by soil formation (paleosols; Fig. 2). Similarly, the thickest unit in the TB succession at PR (90 cm), the fourth unit (CF131/132) above the CI, chemically corresponds to the distal Y-3 tephra (Fig. 2; Fig. S2).

The features of the S-19 drill core (25.9-22.50 m) and PR (CF131-132) units are similar. Both are characterised by a lower sub-unit of fall containing a bed of moderately sorted pumice lapilli (5-10 cm thick), and an upper sub-unit (up to 3.30 m thick in S-19) of poorlysorted ash-dominated pyroclastic density current (PDC) deposits. The PDC sub-unit contains abundant, distinctively large (up to 3 cm diameter) accretionary lapilli which are exclusive to these TB units. The main difference between the two localities is the thickness; the PDC sub-unit of S-19 is considerably thicker than at PR (50cm), either reflecting paleotopography or post-depositional erosion. Detailed descriptions of the S-19/PR tephra successions are provided in Data Repository 1.

The compositional glass data for the PR (samples CF131 and CF132) and S-19 (samples 25.8m and 22.9m) units verify their lithostratigraphic correlation. They are the only TB eruption deposits to contain three compositional groupings (Fig. 2A). The Component 1 glasses are high-SiO 2 (~61.6-62.4 wt.%) trachytic glasses with ~8.3-8.7 wt.% K 2 O (Fig. 2A), they are enriched in incompatible trace elements (e.g., ~30-33 ppm Th), and match the trachytic end-member of the distal Y-3 (Fig. 2). Component 2 glasses are lower in SiO 2 (~59.8-60.6 wt.%), have higher K 2 O contents (~9.7-10.3 wt.%; Fig. 2), and lower levels incompatible trace element enrichment (e.g., ~14-15 ppm Th) -consistent with the distinctive phono-trachytic end-member of the Y-3 tephra (Fig. 2).

The third, less abundant and more chemical variable component, are lower SiO 2 trachytic glasses (~57.6-60.8 wt.%) and show the lowest K 2 O contents (~7.9-8.3 wt.%; Fig. 2) observed in the investigated deposits.. The Component 3 glasses are not observed in the upper portion of the thicker S-19 borehole (25.8-22.50 m) unit, or the widespread Y-3 tephra, this implies the widespread ash dispersal is not associated with the earliest phase of the eruption.

40 Ar/ 39 Ar GEOCHRONOLOGY Sanidine crystals from three near-source pyroclastic samples stratigraphically and chemically related to the distal Y-3 were dated using 40 Ar/ 39 Ar methods at the Laboratoire des Science du Climat et de l'Environment (LSCE), France (sample 22.9m -top of the S-19 borehole unit) and at the Berkeley Geochronology Centre (BGC), USA (sample 25.8mbase of the S-19 unit; sample CF132 upper portion of the PR unit). Weighted mean ages are calibrated to the age of the Alder Creek Sanidine (ACs) =1.1891 ± 0.0008 Ma (1σ, [START_REF] Niespolo | Intercalibration and Age of the Alder Creek 40 Ar/ 39 Ar standard[END_REF] and calculated using the decay constants of [START_REF] Renne | Joint determination of 40 K decay constants and 40 Ar*/ 40 K for the Fish Canyon sanidine standard, and improved accuracy for 40 Ar/ 39 Ar geochronology[END_REF]. Analytical methods and instrumentation used are fully outlined in the Data Repository 1, along with probability diagrams (Fig. S3-7), and R-values (Table S3). The 40 Ar/ 39 Ar weighted mean ages are presented in a summary age probability diagram at 2σ uncertainty (Fig. 3).

The ages of individual crystals analysed from the S-19 22.9 m sample show a bimodal distribution with a dominant younger population (~29 ka) and an older population at ~40 ka (Fig. 3). The 40 Ar/ 39 Ar weighted mean age determined for this upper portion of the eruption deposit is 29.0 ± 0.8 ka [2σ] (MSWD=0.4, P=1, n=18). The older crystals cluster about the age of the CI eruption (40 ka) and are considered xenocrysts from parts of the crystallised magmatic system or deposits of the CI. The whole 2 sigma age range of the CI xenocrysts identified in the S-19 22.9 m sample are used, in conjunction with the stratigraphic position of all correlated samples above the CI, as a basis to exclude CI xenocrysts from the remaining analyses.

The exclusion of CI xenocrysts produces a 40 Ar/ 39 Ar weighted mean age of 28.4 ± 2.1 (MSWD=0.52, P=1, n=44) for the 25.8m basal sample of the borehole unit. The combined 40 Ar/ 39 Ar weighted mean age for the S-19 25.9-22.5m (samples 25.8m and 22.9m) eruption deposit is 29.0 ± 0.8 ka. The CI xenocrysts excluded 40 Ar/ 39 Ar age determinations for the PR deposit (sample CF132), are 31.6 ± 2.8 (MSWD=0.33, P=1, n=23) using Nexus, and 31.3 ± 2.6 (MSWD=0.15, P=1, n=44) from Noblesse, these ages are statistically indistinguishable from one another. The combined 40 Ar/ 39 Ar weighted mean age for the PR eruption deposit (CF131/132) is 31.4 ± 1.9 ka.

The S-19 and PR 40 Ar/ 39 Ar ages are statistically indistinguishable at the 95.4% confidence level and can be combined to provide a 40 Ar/ 39 Ar weighted mean age of 29.3 ± 0.7 ka (MSWD=0.5, P=1, n=96) for the eruption deposit, which perfectly overlap with radiocarbon age determinations for the distal Y-3 tephra (Fig. 3), verifying the stratigraphic and chemical correlations.

ASSESSING THE ASH DISPERSAL AND PDC VOLUMES

The ash dispersal was simulated using a semi-analytical dispersal model HAZMAP [START_REF] Macedonio | A computer model for volcanic ash fallout and assessment of subsequent hazard[END_REF]. Integrating near-source and distal (Y-3) tephra fall deposit thicknesses (Table S1) within a HAZMAP model allows us to constrain the eruption parameters (Table S4; [START_REF] Costa | Tephra fallout hazard assessment at the Campi Flegrei caldera (Italy)[END_REF]. Some reported Y-3 localities could not be used because either layer thicknesses are not reported, or they are likely inflated due to the depositional setting (e.g. marine canyons).

The modelled dispersal pattern for the Y-3 tephra (Fig. 1A) reproduced the recorded thicknesses within an order of magnitude, with all but one of the modelled thickness values being between 1/3 and 3 times the observed thickness at each locality (Fig. S8). The limitations of this modelling approach are outlined in [START_REF] Matthews | Ultradistal tephra deposits from super-eruptions: Examples from Toba, Indonesia and Taupo Volcanic Zone[END_REF], and the uncertainties are large due to the lack of data and assumptions (e.g. deposit density does not change over the area, assumed bulk granulometry). However, these results place constraints on the order of magnitude of the eruption, revealing that the fall comprised ~16 km 3 DRE (assuming a rock density of 2350 kg/m 3 ) of ejected magma, from a column that reached ~59 km high. The volume of the near-source PDC was estimated at a minimum of ~1 km 3 DRE using the Delaunay Triangulation method of [START_REF] Macedonio | An algorithm for the triangulation of arbitrarily distributed points: Applications to volume estimates and terrain fitting[END_REF], which is suitable for reconstructing the volume of geological horizons where the thickness values are irregularly-spaced (Table S1). Combining the HAZMAP fall volume with that of the PDC indicates that the total volume of the eruption would be ~17 km 3 DRE, which equates to a VEI 6 (after [START_REF] Newhall | The Volcanic Explosivity Index (VEI) e an estimate of the explosive magnitude for historical volcanism[END_REF] or at least a magnitude 6 eruption (M6.6, following [START_REF] Pyle | Sizes of volcanic eruptions[END_REF].

THE MASSERIA DEL MONTE TUFF ERUPTION DYNAMICS AND HAZARDS AT CFc

The modelled eruption parameters (Table S4) would imply that the herein named Masseria del Monte Tuff (named after the most proximal TB succession), responsible for the Y-3 tephra, is similar in magnitude and more widely dispersed as a visible fall deposit than the younger NYT (Fig. 1A) caldera-forming eruption (recalculated to ~25 km 3 DRE; following [START_REF] Deino | The age of the Neapolitan Yellow Tuff caldera-forming eruption (Campi Flegrei caldera -Italy) assessed by 40 Ar/ 39 Ar dating method[END_REF] of CFc. Given the amount of magma erupted during the Masseria del Monte Tuff it is likely to have generated a caldera. The lack of physical evidence implies that if one was generated, it must have been within the collapse structure of the younger NYT.

Although the Masseria del Monte Tuff and NYT eruptions were a similar order of magnitude, there is a stark contrast between their PDC deposits. At similar distances from the caldera wall, the PDC deposits of the eruption are thin (~ 3.3 m) relative to the ~16 m thick deposits of the NYT. Consequently, it seems likely that the eruption mechanisms of these two large magnitude events were very different; however, due to the limited near-source exposures of the Masseria del Monte Tuff we must use the distal ash fall and ash dispersal modelling to gain a better understanding of the eruption mechanisms.

Given the estimated magnitude of the eruption, the area covered by the Y-3 tephra is somewhat restricted, illustrated by its absence in the eastern Mediterranean marine records (Fig. 1), this suggests efficient aggregation of fine ash. This premature fallout of the finest material, coupled with the high density of the aggregates (modelled to be 1000 kg/m 3 ) and the presence of large accretionary lapilli in proximal outcrops is consistent with the presence of liquid water in the plume [START_REF] Mastin | Adjusting particle-size distributions to account for aggregation in tephra-deposit model forecasts[END_REF]. The modelled column height and the effective diffusion (accounting for the effects of both atmosphere and volcanic plume turbulence) of ~7x10 4 m 2 /s both imply a large mass flow rate (Costa et al., 2013;[START_REF] Costa | Understanding the plume dynamics of explosive super-eruptions[END_REF], which is consistent with a co-ignimbrite plume in the incipient collapsing regime (estimated between 10 9 and 10 10 kg/s by [START_REF] Costa | Understanding the plume dynamics of explosive super-eruptions[END_REF]. These plumes are transitional between a Plinian and co-ignimbrite column and can extend up to 60 km (cf. [START_REF] Costa | Understanding the plume dynamics of explosive super-eruptions[END_REF]. The modelled results are consistent with the lack of a prominent coarse-grain deposit close to the vent, which suggest an unstable eruption column intermittently generating PDCs. The absence of one of the glass populations in the upper portion of the near-source deposits and in the distal Y-3 implies that the widely dispersed ash was largely associated with a coignimbrite plume generated by late PDC phase/s. Similarly, the co-ignimbrite plume significantly contributed to the enormous CI ash dispersal across the eastern Mediterranean [START_REF] Smith | Tephra dispersal during the Campanian Ignimbrite (Italy) eruption: implications for ultra-distal ash transport during the large caldera-forming eruption[END_REF]. From a hazard perspective, recognition of a large magnitude eruption between the CI (40 ka) and NYT (15 ka) eruptions, drastically reduces the reoccurrence interval of large magnitude events at CFc. With the Masseria del Monte eruption dated at 29 ka, the time elapsed between large magnitude eruptions (≥M6/VEI 6) at this densely populated volcano is halved to ~12.5 kyr.

CONCLUSIONS

The Masseria del Monte Tuff erupted at CFc is dated at 29.3 ± 0.7 ka [2σ] and was a large magnitude (M6.6 or VEI 6) event that erupted around ~17 km 3 (DRE) of magma based largely on the distal ash fall covering an area >150,000 km 2 . Reconstructed eruption parameters using the HAZMAP model suggest the deposits are consistent with dispersal from a co-ignimbrite plume in the incipient collapsing regime. The recognition of such a large-magnitude event, only 29 ka ago, illustrates that the near-source record of volcanism at one of the most productive volcanoes in Europe is not fully representative. It also highlights that by integrating distal and near-source records and using ash dispersal models it is possible to better constrain the magnitude, eruptive mechanisms and tempo of volcanism. This is particularly pertinent at volcanoes that commonly experience caldera collapse and the deposition of thick pyroclastic successions, which make the record of past activity hard to access.

Figure 1 :

 1 Figure 1: (A) The distribution of Y-3 tephra occurrences across the central Mediterranean (TableS1), and for comparison the Neapolitan Yellow Tuff(NYT). Overlain is the HAZMAP isopach map of the Masseria del Monte Tuff linked to the Y-3 tephra. (B) A map of Campi Flegrei caldera (CFc) and it major structures, and sampling localities are shown.

Figure 2 :

 2 Figure 2: (A)The distal Y-3 tephra has two dominant components, a trachytic (Component 1) and phono-trachytic (Component 2) end-member, both consistent with CFc glasses[START_REF] Smith | Tephrostratigraphy and glass compositions of post-15 kyr Campi Flegrei eruptions: implications for eruption history and chronostratigraphic markers[END_REF][START_REF] Tomlinson | Geochemistry of the Phlegraean Fields (Italy) proximal sources for major Mediterranean tephras: implications for the dispersal of Plinian and co-ignimbritic components of explosive eruptions: Geochim[END_REF]. The TB units related here to the distal Y-3 at Ponti Rossi (PR) (CF131/132) and in the S-19 Borehole (25.9-22.5m) have the same two components, plus a population of lower SiO 2 trachytes (Component 3). (B) The two endmembers of the distal Y-3 are distinguished using their K 2 O content of their volcanic glasses (grey boxes), only the third unit in borehole S-19, and the third (CF129) and forth unit (CF131-132) in the PR sequence have both component 1 and 2 compositions. The thin ash unit (3; CF129) underlying Unit 4 at PR is inconsistent with the Y-3 tephra based on its levels of trace element enrichment (Fig.S1E-F). (C) Trace elements analysis verifies that both end-members of the Y-3 (Component 1 and 2) are observed in the thickest near-source units at PR and in the S-19 borehole.

Figure 3 :

 3 Figure 3: 40 Ar/ 39 Ar ages (2σ) of proximal CFc samples of the Masseria del Monte Tuff, and terrestrial (Albert et al., 2015) and marine (Ramsey et al., 2015) 14 C ages of the distal Y-3. Individual sanidine age determinations are shown with 1σ errors and xenocrysts (mostly 40 ka CI crystals) are removed from all age calculations.

ACKNOWLEDGEMENTS

PGA and VCS acknowledge funding from Leverhulme (ECF-2014-438) and the John Fell Fund. PRR and EMN thank the Ann and Gordon Getty Foundation and Chihiro Ishida for sample preparation assistance. AC and RI acknowledge the MIUR project Premiale Ash-RESILIENCE Dr. Emma Tomlinson and John Pouncett are thanked for assistance with the trace element analysis, and cartography, respectively. We thank Martine Paterne, Jenni Barclay and a anonymous Reviewer for their constructive comments.