Alexandre Carrara
email: carrara.alexandre.univ@gmail.com

Alain Burgisser

George W Bergantz

The architecture of intrusions in magmatic mush

Keywords: 3 Tables, 50 references Mush, Magma, Intrusion, Density contrast, CDF-DEM, Granular mechanics

published or not. The documents may come

Introduction:

Evidence for injections of new magmas, also called recharge events, are ubiquitous in magmatic systems [START_REF] Wiebe | Mafic replenishments into floored silicic magma chambers[END_REF]. They are inferred to cause the formation of long-lived, supersolidus magmatic reservoirs located in the upper crust (e.g. [START_REF] Annen | Construction and evolution of igneous bodies: Towards an integrated perspective of crustal magmatism[END_REF][START_REF] Annen | The Genesis of Intermediate and Silicic Magmas in Deep Crustal Hot Zones[END_REF][START_REF] Dufek | Lower Crustal Magma Genesis and Preservation: a Stochastic Framework for the Evaluation of Basalt-Crust Interaction[END_REF][START_REF] Karakas | Lifetime and size of shallow magma bodies controlled by crustal-scale magmatism[END_REF]. Together with the thermal structure of the upper crust and the frequency of recharge, the geometry and mode of emplacement of the intruded magma was also identified as having a crucial effect on the long-term evolution of igneous bodies [START_REF] Annen | Construction and evolution of igneous bodies: Towards an integrated perspective of crustal magmatism[END_REF]. Diverse evidence supports the view that magmatic reservoirs reside most time in a mush state that is frequently disturbed by injection of new magmas (e.g. [START_REF] Bachmann | Silicic magma reservoirs in the Earth's crust[END_REF]Cashman et al., 2017, and references therein). A magmatic mush is a crystal-rich magma in which crystals are in close and sometimes frictional contacts, forming a semi-rigid framework where stress is transmitted by force chains [START_REF] Bergantz | On the kinematics and dynamics of crystal-rich systems[END_REF]. As a result, mushes transition between crystal-rich suspensions to a 'lock-up' state that inhibits the ability of the magma to erupt.

The injection of hotter magma into a cooler host has been suggested as a means to trigger volcanic eruptions (e.g. [START_REF] Caricchi | Frequency and magnitude of volcanic eruptions controlled by magma injection and buoyancy[END_REF] and the intrusion style plays a fundamental role in the way mush rejuvenates (process of recycling the mush to generate an eruptible magma) prior to eruption (Parmigiani et al., 2014, and references therein). Several scenarios assume that the intruder is emplaced as sills at the base of the mush (underplating), and rejuvenate it by supplying heat but no mass except possibly exsolved volatiles [START_REF] Bachmann | Gas percolation in upper-crustal silicic crystal mushes as a mechanism for upward heat advection and rejuvenation of near-solidus magma bodies[END_REF][START_REF] Bergantz | Underplating and Partial Melting: Implications for Melt Generation and Extraction[END_REF][START_REF] Burgisser | A rapid mechanism to remobilize and homogenize highly crystalline magma bodies[END_REF][START_REF] Couch | Mineral disequilibrium in lavas explained by convective self-mixing in open magma chambers[END_REF][START_REF] Huber | Thermo-mechanical reactivation of locked crystal mushes: Melting-induced internal fracturing and assimilation processes in magmas[END_REF]. Other scenarios consider that the injected magma may penetrate the mush, producing various degrees of mixing with the resident mush depending on its buoyant acceleration (e.g. [START_REF] Bergantz | Non-stationary entrainment and tunneling eruptions: A dynamic link between eruption processes and magma mixing[END_REF][START_REF] Koyaguchi | Thermal evolution of silicic magma chambers after basalt replenishments[END_REF][START_REF] Weinberg | Mingling in mafic magma chambers replenished by light felsic inputs: fluid dynamical experiments[END_REF]. Whether an intrusion generates extensive mass transfer, or is limited to thermal exchanges between an underplated intruder and a host mush is thus a key element shaping the outcome of open-system events. A major obstacle to our current understanding of the formation and evolution of igneous bodies is that little is known about the architecture of intrusions and controlling physical parameters.

Traditionally, mush rejuvenation scenarios have been based on the results of experiments performed with pure fluids mimicking the bulk physical properties (density and viscosity) of the magmas (e.g. [START_REF] Huppert | Replenishment of magma chambers by light inputs[END_REF][START_REF] Jellinek | Mixing and compositional stratification produced by natural convection: 2. Applications to the differentiation of basaltic and silicic magma chambers and komatiite lava flows[END_REF][START_REF] Snyder | Replenishment of magma chambers: comparison of fluid-mechanic experiments with field relations[END_REF]. Mush dynamics, however, differs from that of pure fluids because of the complex rheological feedbacks between melt and crystals. An essential physical process is that melt and crystals may experience relative motions. Numerical simulations explicitly accounting for such decoupled motions as well as the building and destruction of force chains between crystals [START_REF] Bergantz | Open-system dynamics and mixing in magma mushes[END_REF][START_REF] Schleicher | Time scales of crystal mixing in magma mushes[END_REF][START_REF] Schleicher | The Mechanics and Temporal Evolution of an Open-system Magmatic Intrusion into a Crystal-rich Magma[END_REF] have revealed that the local injection of pure melt of the same density and viscosity as the mush interstitial melt easily fluidizes, penetrates, and partially mixes with the overlying mush if it is sufficiently vigorous. This local unlocking of a mush shows that the conditions for efficient mass transfer and mixing are easier to achieve than previously thought. Conversely, it is adding constraints on rejuvenation scenarios based on the emplacement of an underlying mafic gravity current suggesting that underplating may require contrasts in densities and/or viscosities to hinder fluidization.

Our capacity to interpret the various natural expressions of open-system events, such as eruptive products containing both the intruded magma and the resident mush, is hindered by our partial understanding of the architectural end-members of these events, such as fluidization or underplating. To characterize the geometry and emplacement styles of intrusion events into a residing mush, we performed numerical simulations using a combination of fluid mechanics and discrete elements [START_REF] Bergantz | Open-system dynamics and mixing in magma mushes[END_REF][START_REF] Schleicher | Time scales of crystal mixing in magma mushes[END_REF][START_REF] Schleicher | The Mechanics and Temporal Evolution of an Open-system Magmatic Intrusion into a Crystal-rich Magma[END_REF];[START_REF] Mcintire | On the hydrodynamics of crystal clustering[END_REF]Burgisser et al., 2020) . As the dissimilarities between the density and viscosity of the two melts require special attention to better characterize the end-member cases of open-system events, we explored how these parameters condition the dynamics of the intruded material when injected into a mush. We first introduce the numerical model and the dimensionless parameters controlling recharge dynamics that are varied in the simulations.

Results of numerical simulations involving magmas of contrasting physical properties are then presented in the framework of the dimensionless parameters. Finally, we relate our results to well-documented cases of eruptions triggered by an intrusion event.

2: Method

In order to characterize the geometry and emplacement mechanism of intrusion in mush accounting for granular dynamics, we performed Computational-Fluid-Dynamic and Discrete-Element-Method (CFD-DEM) numerical simulations by using the MFIX-DEM software (https://mfix.netl.doe.gov/). The coupled melt-crystal flow is computing solving the mass and momentum conservation with a finite volume method for the fluid phase and a discrete model including explicitly of particle's interactions for the crystal phase. Details about the theory and implementation of the model can be found in [START_REF] Garg | Open-source MFIX-DEM software for gas-solids flows: Part I-Verification studies[END_REF], [START_REF] Syamlal | MFIX documentation numerical technique[END_REF], [START_REF] Syamlal | MFIX documentation theory guide[END_REF], and validation of the DEM approaches in [START_REF] Garg | Open-source MFIX-DEM software for gas-solids flows: Part I-Verification studies[END_REF] and [START_REF] Li | Open-source MFIX-DEM software for gas-solids flows: Part II-Validation studies[END_REF] (see supplementary information 1 for an overview of the numerical model, a discussion about its validity, and a list of the equations we used).

The computational domain is a 3D medium of 1.6 × 0.8 × 0.05 m (length × height × width; 128 × 64 × 4 cells) filled with a resident mush (Fig. 1). This geometry also allowed us to populate the mush with mm-size particles, ensuring that the dimensionless numbers describing the nature of the flow and the contrasts between the physical properties of the two magmas cover the ranges found in nature (see Supplementary information 1). We will show a posteriori that our particle bed behaves identically to a bed twice as thick [START_REF] Bergantz | Open-system dynamics and mixing in magma mushes[END_REF]. Our runs are thus representative of an open system event despite the small size of the domain compared to a natural system. We used such geometry instead of a two dimensional one to ensure that the build-up and breaking of force chains have a sufficient degree of freedom in space to replicate best the mechanics of the granular phase. We created a mush layer of ~0.3 m height with an initial crystal volume fraction of ~0.64 by simulating the settling of the particles in a vacuum and positioning them at the base of the domain. We used the same density for all particles ( kg m -3 ) and three different diameters (4.5, 5, and 5.5 mm, constituting 25, 50, and 25 % of the total number of particles, respectively) to avoid artificial self-organization of the particles in a regular lattice. All simulations use the same initial particle bed. A crystal-free magma is injected at the base of the mush layer with a superficial vertical velocity, , through an inlet having a width, . The density and the viscosity of the injected melt are kept constant between all the simulations ( kg m -3 ; Pa s, see table 2 for the list of the parameters kept constant). We used a conduit of 3.2 cm in height to supply the inlet to ensure that the intruder enters the mush as a Poiseuille flow. At the top of the domain, we used a pressure outflow boundary condition to ensure the overall mass conservation within the entire domain, which is consistent with an opensystem event. The boundary conditions at the front and back of the domain are cyclical, which means that the intruder corresponds to a dyke having one infinite dimension. All the other boundary conditions are non-slip walls (Fig 1). To maintain constant values of melt density and viscosity during the runs (and hence constant density and viscosity contrasts), thermal effects are ignored. This is consistent with the small dimensions of the computational domain that ensure run times shorter than those allowing significant heat exchanges (the characteristic time to diffuse the heat over 1 cell is ~400s and the longest simulations lasts ~500s) In each fluid cell, the local melt density and viscosity depend linearly on the two end-member properties and on their respective concentrations, which are tracked by a transport equation (see Supplementary material 1).

We performed simulations by varying the density and viscosity of the host melt. In order to compare simulations, we used dimensionless quantities to scale the effects of the contrasts in densities and viscosities, and injection velocities. The injection velocity and melt viscosity control the stress applied by the input of new materials to the mush. These parameters enter the minimum fluidization velocity, [START_REF] Schleicher | Time scales of crystal mixing in magma mushes[END_REF], see supplementary information 2 for derivation of ), which expresses the superficial velocity at which the upward drag force equals the reduced weight of the particle bed. This threshold is proportional to the density contrast between the particles and the surrounding melt and inversely proportional to melt viscosity. As the injected melt differs from the host melt, two minimum fluidization velocities can be calculated depending on which melt is considered. For all simulations, we used the minimum of these two velocities, which here always corresponds to that using the host melt properties. The dimensionless injection velocity, , is defined as:

.

In simulations having identical , the injection imposes the same stress to the overlying mush.

However, the time needed to inject the same new melt volume changes between simulations because varies. We thus used a dimensionless time, , to scale the simulation time [START_REF] Bergantz | On the kinematics and dynamics of crystal-rich systems[END_REF]:

, ( 2
)
where is the simulation time. In this way, simulations having identical implies that the same volumes of intruder have been injected until that dimensionless time and simulation results can be compared directly. We use the Atwood number to scale the density contrast between the two materials. Here, two Atwood numbers may be defined. The first one, , expresses the buoyancy contrast between the two melts:

, ( 3
)
where is the density of the intruded melt, and is the host melt density. The second one, , takes the presence of crystals in the host material into account and scales the bulk densities (the weight of material including crystals and melt per unit of volume):

, ( 4
)
where is the density of the host solids, and is the particle volume fraction. A negative Atwood number indicates that the intruder is buoyant compared to the mush, whereas a positive one indicates a tendency to sink. The viscosity contrast, , between the two melts is expressed as:

, ( 5
)
where is the host dynamic viscosity and is that of the injected melt.

3: Results

We performed 25 numerical simulations to explore the influence of the host melt density and viscosity (See Table 3 for a list of all the simulations and corresponding parameters). For these simulations, the injection velocities are such that the ratio with the respective minimum fluidization velocity, , remains constant at =21.2. This ratio is chosen to match that used previously in similar works [START_REF] Schleicher | Time scales of crystal mixing in magma mushes[END_REF][START_REF] Schleicher | The Mechanics and Temporal Evolution of an Open-system Magmatic Intrusion into a Crystal-rich Magma[END_REF] according to the formula presented in the supplementary material 2. We performed an additional 4 simulations at higher injection velocities to explore the effect of on intrusion dynamics. The fluidization regime was observed in the simulations once , and consists in the development of a fluidized area above the inlet in which the intruded melt rises through the mush (Fig. 3A-C), as described previously [START_REF] Bergantz | Open-system dynamics and mixing in magma mushes[END_REF][START_REF] Schleicher | Time scales of crystal mixing in magma mushes[END_REF]. The fluidization of the mush is initiated by the dilation of the crystal framework to crystal volume fraction below 0.3 above the inlet that locally destabilizes the forces chains network that supports the bed and separates the crystals in contact. The fluidized volume grows vertically above the inlet because of two mechanisms. The first is the upward entrainment of the particles localized above the fluidized cavity, which results in bulging the top surface of the mush layer (Fig. 3A-C). The second mechanism is the progressive erosion of the crystals jammed at the boundary between the mush and the fluidized volume. Once separated, crystals start settling in the fluidized area because of this process of mush erosion, causing the fluidized area to ascend faster than the intruded melt (green outline in Fig 3A -C). The intruder flows mainly vertically with a minor lateral porous flow. When the fluidized cavity reaches the top of the particle bed, its width progressively decreases to stabilize in the shape of a vertical chimney. At steady state, when t * >1, the crystals located within the chimney show both upward and downward motions whereas the ones located around the chimney flow slowly in the direction of the inlet, forming a 'mixing bowl' as a whole, fully recovering the dynamics first described in [START_REF] Bergantz | Open-system dynamics and mixing in magma mushes[END_REF].

The spreading regime, which prevails in simulations once , is characterized by the lateral spreading of the injected melt similarly to a gravity current hugging the floor of the host reservoir (Fig. 3D-F). The main difference with a pure fluid gravity current is that the melt is progressively flowing across the mush as permeable flow. At the start of the injection, the crystal framework experienced a dilation, which initiates host crystals settling in the same fashion as in the fluidization regime. The lateral flow of the intruded melt is able to laterally entrain the host crystals, creating two counter rotating granular vortexes in the residing mush with downward motions above the inlet (Fig. 3D-F). Such granular vorticity affects the flow pattern of the fluid in the mush. The fluidized volume grows either predominantly laterally or vertically, depending on the relative importance between the lateral entrainment of the host solids by the intruder and the vertical settling of the mush crystals. As the lateral propagation of the intruder progresses, so does the size of the two granular vortexes, making this style of intrusion affect a larger mush volume than the fluidization regime.

The rising regime (Fig. 3G-I), is characterized by the ascent of the intruded melt within the mush that occurred in simulations once . Runs start with the initial growth above the inlet of a cavity filled with the intruded fluid. The cavity becomes gravitationally unstable and ascends within the mush, forming a Rayleigh-Taylor instability. The ascent of the intruder continues above the particle bed, entraining solids from the host. The instability is driven by its head because of the buoyant batch of intruded melt. This driving batch is surrounded by a volume of fluidized host mush (Fig. 3G-I, better seen in Supplementary Movie 3 as the absence of particle overlap). The dimensionless time at which the intrusion reaches the mush top (t * ~0.3) is shorter than that of the two other regimes because the Rayleigh-Taylor instability significantly accelerates the transport of the intruder.

Figure 2 suggests that the viscosity contrast does not control the end-member shape of the intruder flow. Larger viscosity contrasts, however, increase the trends of some aspects of mush dynamics. Figure 3 illustrates how viscosity bears on flow patterns.

In the fluidization regime, the increase of the host viscosity enhances the formation of crystal-poor batches at the top of the intruded volume (Fig. 3A-C). Because the minimum fluidization velocity within the intruded melt is higher than for the host, the crystals are not fluidized and sediment in the intruded melt to accumulate atop the inlet (Fig 3B -C). Because we defined t * to scale the dynamics of the mush, the increase of the host melt viscosity decreases the injection velocity and increases the duration, , required to reach the dimensionless time t * =1. As a result, increasing melt viscosity decreases the ability for the intruded melt to experience lateral porous flow through the host crystal frameworks (Fig. 3B-C). It also increases the ratio between crystal terminal velocity in the intruded and host melts, which results in the formation of the crystal poor volume at the top of the intruded volume (Fig. 3B-C). The increase in the host melt viscosity, however, does not affect the volume of mush inflated by the injection and showing distortion of the force chains.

In the spreading regime, high viscosity contrasts enhance the lateral spreading of the intruder and the entrainment of the host crystals in the two counter rotating vortexes (Fig 3E -F).

Large host melt viscosity causes the lateral entrainment of the solids to be more efficient than does not significantly vary with the host melt viscosity (Fig. 3G-I).

The additional 4 simulations in the spreading regime suggest that buoyancy effects dominate the intruder flow up to . Figure 4 shows the temporal evolution of the height reached by the intruded volume, , as a function of injection rate. All injections grow purely vertically at first ( ). As seen above, at the low injection rate of 21.2, the intrusion stalls rapidly and spreads laterally (simulation A25, Fig. 4). Increasing the injection rate causes stalling to occur later and higher. When , injection growth switches from vertical to radial. When , the behavior of the intruder is dominated by the injection rate, which causes the radially growing intrusion to reach the top of the bed at . Despite that all simulations have the same intruder shape before stalling, the size of the region surrounding the intruder that is affected by dilatancy increases with . The highest injection rate (simulation B4 with ) strictly follows the theoretical curve for a radial growth and reaches at , as predicted by geometrical arguments (supplementary information 3).

Mush dilation is key for permeable melt flow to occur. The initiation of the intrusion increases the pore pressure in the mush around the inlet (Fig. 5A). This overpressure progressively propagates outwards and decreases the crystal volume fraction in the overlying mush (Fig. 5B). As the intrusion propagates, the effect of the overpressure is supplemented with the Reynolds dilatancy generated by the granular vortexes in the mush (Fig 5C). The dilation of the solid framework increases its permeability of the solid framework and in turn the possibility of relative motion between the crystals and the interstitial melt (Fig. 5C). This phenomenon is particularly clear in the case of the rising regime. The intruder is surrounded by a volume of mush that underwent such dilation that it is in the dilute regime. The contact region between the two magmas is dominated by melt-melt interface interspersed with isolated crystals. As a result, entrainment is ruled by melt vorticity.

4: Comparison with natural systems

To test the applicability of our results to natural cases, we gathered from the literature the physical parameters of 13 eruptions involving the intrusion of new magma (Table S3-S4 in the supplementary information 5). All host magmas are well-characterized but for a few cases for which there is ambiguity on the respective roles of the intruder and host magmas (Unzen, Minoan, and Katmai-Novarupta). In the studies surveyed, melt viscosity and melt density of host magmas were most often directly determined from eruptive products and pre-eruptive conditions such as pressure, temperature, and melt water content (details on how parameters were obtained are in Table S3-S4 (see supplementary information 5)).

Cases are organized into three categories depending on the observed eruptive sequence. In the first category, the intruder was erupted first, followed by the emission of host magma or a mixture of host and intruder. This category implies that the intruder magma was able to efficiently penetrate and pass through the host magma. In the second category, both host and intruder magmas were erupted simultaneously, with the intruder most often forming enclaves or mingling structures. The last category feature cases where the mixing was so thorough that the eruptive products only bear cryptic traces of the intruder, such as isolated intruder crystals floating in the host or crystal disequilibrium textures.

Figure 6A shows the ratios of bulk viscosities and bulk densities between the intruder magma and the host magma(s) for the 13 eruptions. Figure 6A contains two physically meaningful thresholds, that of neutral buoyancy at the bulk density ratio of 0 and that of equal viscosity at the bulk viscosity ratio of one. The three types of eruptive sequence are not sorted following any of these thresholds. Figure 6B shows the same eruptions plotted as functions of melt properties instead of bulk properties. Our numerical runs cover the full range of natural density ratios and a more restricted range of viscosity ratios (from 1 to 10 2 vs. 10 -2 to 10 4 in nature). Figure 6B also shows the dividing line between rising and spreading dynamics at the level of neutral buoyancy with respect to the melts. With the possible exception of the 1912 Katmai-Novarupta eruption (see Discussion), the rising regime is populated by the eruptions that first ejected intruder material. This divide between cases where at least some of the intruder magma had the capacity to go unscathed through the host and cases where none of it escaped from host interaction is consistent with our numerical results.

5: Discussion

Our results are helpful to predict the behavior of an intrusion within a mush. The Atwood number between the two melts, , is the parameter having a first order control on the geometry of the intrusion. On the contrary, the commonly used level of neutral bulk buoyancy (e.g. [START_REF] Huppert | Replenishment of magma chambers by light inputs[END_REF][START_REF] Snyder | Replenishment of magma chambers: comparison of fluid-mechanic experiments with field relations[END_REF] does not mark any particular change in dynamic behavior (Fig. 2). This result illustrates that the relative motion between the solids and surrounding melt is of primary importance when studying mush processes. Experiments, or numerical simulations, that account of the presence of the solids or exsolved volatiles as discrete entities (e.g. [START_REF] Barth | Valve-like dynamics of gas flow through a packed crystal mush and cyclic strombolian explosions[END_REF][START_REF] Bergantz | Open-system dynamics and mixing in magma mushes[END_REF][START_REF] Girard | Buoyant replenishment in silicic magma reservoirs: Experimental approach and implications for magma dynamics, crystal mush remobilization, and eruption[END_REF][START_REF] Hodge | Experimental constraints on the deformation and breakup of injected magma[END_REF][START_REF] Mcintire | On the hydrodynamics of crystal clustering[END_REF][START_REF] Michioka | Rayleigh-Taylor instability of a particle packed viscous fluid: Implications for a solidifying magma[END_REF][START_REF] Parmigiani | Mush microphysics and the reactivation of crystal-rich magma reservoirs[END_REF][START_REF] Schleicher | Time scales of crystal mixing in magma mushes[END_REF][START_REF] Schleicher | The Mechanics and Temporal Evolution of an Open-system Magmatic Intrusion into a Crystal-rich Magma[END_REF] are the most likely to faithfully reproduce mush dynamics.

Neglecting phase decoupling by considering the magma as a single-phase fluid having effective properties such as bulk density or bulk viscosity will not capture the blending of crystal contents between host and intruder and the simultaneous but independent evolution of the melt-melt interface (Fig. 5C).

The decoupling between the motions of the two phases results from processes unique to granular mechanics that our discrete numerical model is able to capture. Mush dilation causes entrainment by melt vorticity alone (Fig. 5). Efficient entrainment of two fluids with a viscosity contrast occurs only when the most viscous fluid bears large levels of vorticity [START_REF] Jellinek | Mixing and compositional stratification produced by natural convection: 2. Applications to the differentiation of basaltic and silicic magma chambers and komatiite lava flows[END_REF]. In our runs, the intruder melt viscosity is equal or less than that of the host, and the vorticity is concentrated close to or inside the intrusion (Fig. S4 in the supplementary information 4). This situation yields the weak entrainment observed in the rising regime and the transition from vertical growth to spreading of the intrusion melt as injection velocity decreases (Fig. 4).

The concept of bulk density thus fails to predict the intrusion geometry for two reasons. First, it assumes the absence of relative motion and thus ignores the transfer of crystals from host to intrusion. Second, in cases when sufficient mush dilation occurs, entrainment is controlled by the melt-melt interface and the associated density and viscosity contrasts. The interplay between pore pressure, dilation, melt interface dynamics, and permeable flow controls the transport of mass within our modeled magmatic reservoir.

We characterized the parameter ranges of a series of well-documented cases of eruptions that features magma mixing, focusing on the densities and viscosities of the two end-member magmas involved and on the order of the eruptive sequence. One case, Katmai, straddles two eruptive sequence categories because the intruders may have been transported alongside (as opposed to through) the host magmas. Both scenarios are close to the neutral buoyancy level (Fig. 6B) and each individual scenario is consistent with our regimes. The Katmai eruption first emitted rhyolite. The Katmai scenario corresponding to a rhyolite intruding a more mafic host [START_REF] Eichelberger | Eruption of andesite triggered by dyke injection: contrasting cases at Karymsky Volcano, Kamchatka and Mt Katmai, Alaska[END_REF] is consistent with it being located in the rising regime. In the other scenario [START_REF] Singer | The timing of compositionally-zoned magma reservoirs and mafic 'priming' weeks before the 1912 Novarupta-Katmai rhyolite eruption[END_REF], the rhyolite is part of the host reservoir, which is consistent with that scenario being in the spreading regime. The overall good agreement between the observed eruptive sequences and our numerical results (Fig. 6B) supports the fact that opensystem events are, to first order, controlled by the density contrast between the melt phases of the intrusion and mush. It also suggests that injection momentum was quickly exhausted, letting buoyancy control the unfolding of the event.

One special natural case can be evaluated against our dynamic regimes. The intruder of the 1991-1995 eruption of Unzen volcano could have been either andesitic [START_REF] Holtz | Experimental petrology of the 1991-1995 Unzen dacite, Japan. Part I: phase relations, phase composition and pre-eruptive conditions[END_REF], or basaltic [START_REF] Browne | Magma mingling as indicated by texture and Sr/Ba ratios of plagioclase phenocrysts from Unzen volcano, SW Japan[END_REF]. Regardless of its composition, the intrusion caused thorough mixing and the first magma erupted was the product of this mixing. If the intruder was basaltic, it was buoyant with respect to the felsic host and if it was andesitic, it was denser that the host. As a result, Unzen spans the divide between the rising and spreading regime (Fig. 6B). Considering that the intruder input was large (>30 wt% of the eruptive products; Holtz et al., 2005), and if any credit is given to our inferences, the intruder was more likely to be andesitic than basaltic because this latter composition would have been prone to preserve its integrity while going through the host mush, erupting first.

6: Conclusions

This study highlights the importance of granular mechanics in mush processes, which differ significantly from ones expected with purely fluid models. As expected, our simulations show that when the injection velocity is high ( ), intrusion dynamics is dominated by the injection momentum and the intruded cavities grow radially. When the injection velocity is below this threshold, however, buoyancy controls the behavior of the intruder in an unexpected way.

Bulk buoyancy contrasts appear to play no role in the way the intruder flows. Instead, the density contrast between the host and intruded melts exerts a first-order control on the architecture of an intrusion event in a mush. When the two melt densities are identical, the intruder fluidizes the mush and creates a mixing bowl, as described in [START_REF] Bergantz | Open-system dynamics and mixing in magma mushes[END_REF]. When the intruded melt is lighter than that of the host, it rises through the mush. Mush dilation around the intruder causes the contact region between the two magmas to be dominated by melt-melt interface interspersed with isolated crystals. Entrainment in this rising regime is ruled by the amount of vorticity of the most viscous melt. As in our runs the intruder melt viscosity was equal or lower than that of the host, no entrainment was observed. Intruder melts denser than the host spread laterally partly as permeable flows through the host mush. The lateral spreading of the intruder generates two counter rotating granular vortexes with downward motions above the inlet, which maximizes the volume of the mush entrained by the gravity current. In this spreading regime, the combined effects of the initial pore overpressure at the inlet and the Reynolds dilatancy resulting from the lateral spreading of the intruder are able to fluidize the overlying mush.

We tested whether the first-order effect of melt density contrast was expressed in nature.

We tallied 13 well-documented eruptive sequences, classifying them according to the expected outcomes of the three dynamic regimes we defined. We found overall good agreement between eruption sequences and our model predictions, which suggests that pore pressure, dilatancy, and permeable flow play a fundamental role in the unfolding of open-system events. Granular dynamics and the decoupling of melt and crystals are thus key in shaping reservoir and volcanic processes.

Tables:

Symbol (unit) Definition

At

Atwood number between the two melts At b

Bulk Atwood number

Figure 2

 2 Figure 2 plots the simulations at the lowest , 21.2, as functions of the dimensionless

  Fig 3H-I, and a significant volume of fluid is injected before the instability accelerates. On the

1 Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 11234 Figure 1: Simulations initial condition. [A] The drawing represents the computational domain2

Acknowledgments:

All the numerical simulations presented in this paper were performed using the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA: http://www.ci-ra.org). This work was partially supported by the [START_REF] Bergantz | Open-system dynamics and mixing in magma mushes[END_REF]. See text for details regarding the special case of Unzen.

57


Oil quality and content in olives depend on complex agricultural factors that determine the ripening process of fruits [START_REF] Gracia | Non-destructive assessment of olive fruit ripening by portable near infrared spectroscopy[END_REF] . On the other hand, the oil content (OC) determines the adjustment of key parameters of the milling process [START_REF] Armenta | The Use of Near-Infrared Spectrometry in the Olive Oil Industry[END_REF] . The estimation of the average oil content of incoming olive batches in a mill constitutes the basis of oil extraction control in terms of duration and temperature of malaxation [START_REF] Trapani | A kinetic approach to predict the potential effect of malaxation time temperature conditions on extra virgin olive oil extraction yield[END_REF] , and of the rate of feeding in the decanter pump [START_REF] Guerrini | Changes in Olive Paste Composition During Decanter Feeding and Effects on Oil Yield[END_REF] . All of which ensure the optimization of oil extraction: improvement of yield (litres per kg) and quality of olive oil [START_REF] Guerrini | Changes in Olive Paste Composition During Decanter Feeding and Effects on Oil Yield[END_REF] .

Currently, gravimetric analysis under Soxhlet extraction is the official method to determine the OC in olive batches 5 , however, it is time consuming, and requires sample preparation and dissolvents. More recently, nuclear magnetic resonance (NMR) is being used in quality laboratories to determine OC. Both methods make use of olive paste in such environments.

García, et al. [START_REF] García | Comparative study of various analytical techniques (NIR and NMR spectroscopies, and Soxhlet extraction) for the determination of the fat and moisture content of olives and pomace obtained from Jaen (Spain)[END_REF] have shown that either NMR, or Soxhlet extraction provide comparable estimations in oil concentration for milled olives, and thus both methods are taken for redundant.

NIR spectroscopy is an alternative method for oil content quantification increasingly used in quality laboratories, and industries for routine analysis (commercial equipment Foss OliviaTM; Bruker MPA). It is fast and does not require dried samples, but it is most frequently applied on olive paste [START_REF] García | Comparative study of various analytical techniques (NIR and NMR spectroscopies, and Soxhlet extraction) for the determination of the fat and moisture content of olives and pomace obtained from Jaen (Spain)[END_REF][START_REF] Allouche | Near infrared spectroscopy and artificial neural network to characterise olive fruit and oil online for process optimisation[END_REF] . More recently, innovative implementations of NIR systems for OC quantification have been conducted on intact olives, for breeding programs [START_REF] Leon-Moreno | Usefulness of portable near infrared spectroscopy in olive breeding programs[END_REF][START_REF] De La Rosa | Fruit characteristics and fatty acid composition in advanced olive breeding selections along the ripening period[END_REF] where the selection of individuals becomes the main target, or for the evaluation of fruit entering the mill [START_REF] Fernandez-Espinosa | Combining PLS regression with portable NIR spectroscopy to on-line monitor quality parameters in intact olives for determining optimal harvesting time[END_REF][START_REF] Salguero-Chaparro | Near infrared spectroscopy (NIRS) for on-line determination of quality parameters in intact olives[END_REF][START_REF] Salguero-Chaparro | Feasibility of using NIR spectroscopy to detect herbicide residues in intact olives[END_REF] . In addition, advances in NIRS technology have allowed the evolution from laboratory equipment [START_REF] Kavdir | Prediction of olive quality using FT-NIR spectroscopy in reflectance and transmittance modes[END_REF][START_REF] Hernandez-Sanchez | From NIR spectra to singular wavelengths for the estimation of the oil and water contents in olive fruits[END_REF] to the implementation of in-field portable devices [START_REF] Barreiro | In Iv International Symposium on Applications of Modelling as an Innovative Technology in the Agri-Food-Chain: Model-It[END_REF] as well to on-line spectrometers [START_REF] Giovenzana | Use of visible and near infrared spectroscopy with a view to on-line evaluation of oil content during olive processing[END_REF] , leading to faster and more efficient analysis compared to laboratory NIRs. Still, NIR remains an indirect method that requires a rigorous calibration procedure to be implemented.

Moving from mill batches, to small samples or to individual olives requires the consideration of concepts such as homogeneity/heterogeneity from the point of view of the theory of sampling [START_REF] Esbensena | Theory of sampling (TOS)-fundamental definitions and concepts[END_REF] . Different sampling strategies require different management of heterogeneity.

In some studies, batches constitute the decision units with two alternative procedures: using homogenized samples such as pastes (composite samples), or the use of intact fruits in the case of non-destructive methods [START_REF] Salguero-Chaparro | Near infrared spectroscopy (NIRS) for on-line determination of quality parameters in intact olives[END_REF][START_REF] Salguero-Chaparro | Feasibility of using NIR spectroscopy to detect herbicide residues in intact olives[END_REF][START_REF] Giovenzana | Use of visible and near infrared spectroscopy with a view to on-line evaluation of oil content during olive processing[END_REF] .

NMR and chemical extraction protocols (factory standard) [START_REF] Kavdir | Prediction of olive quality using FT-NIR spectroscopy in reflectance and transmittance modes[END_REF][START_REF] Deblangey | Comparative study of two methods (hexane extraction and NMR) for the determination of oil content in an individual olive fruit[END_REF][START_REF] Cayuela | NIR prediction of fruit moisture, free acidity and oil content in intact olives[END_REF] have also been implemented to determine the oil content of individual olive fruit. Comparisons of the official methods with regard to NMR quantification of the oil content in olives has demonstrated that NMR presents the highest overall efficiency (more sensitivity, good repeatability and higher precision) [START_REF] Hernandez-Sanchez | From NIR spectra to singular wavelengths for the estimation of the oil and water contents in olive fruits[END_REF][START_REF] Deblangey | Comparative study of two methods (hexane extraction and NMR) for the determination of oil content in an individual olive fruit[END_REF] . Thus, NMR is a more direct procedure for the calibration of NIR models.

At the industrial level, computer vision is a widely used technology in the production of table olives to determine the fruit size and detect external damage (H2020-SMEInst-2018-2020-2 Project: Evoolution). The possibility of online scanning of all the olives to be processed in a mill, either by multi or hyperspectral image systems (VIS and NIR) [START_REF] Guzman | Determination of the olive maturity index of intact fruits using image analysis[END_REF] , opens the doors to have a very accurate and real-time information of the flow of the oil that effectively enters the industrial process, even allowing previous segregation of individual fruits in more homogeneous batches.

The purpose of this work was to evaluate different strategies to calibrate NIR models in practical situations in which the decision units are individual olive fruits for NIR, and olive paste for the reference method. This paper addresses, as an innovation compared to previous literature, three approaches managing different levels of heterogeneity in the references considered to estimate oil concentration (%) and total oil content (g).

Materials and Methods

Sampling and methods. Olive fruits from a commercial mill in Toledo (D.O. Montes de Toledo), belonging to the varieties Arbequina, Picual and Cornicabra, were harvested with a wide maturity range at 12 harvest dates from November 2015 to February 2016. In total, 278 olives were clustered in 31 batches according to the harvest date, variety and maturity level (regarding to the external colour of the olives: green, purple and black).

Once the olives were harvested and classified into 31 batches of similar maturity, they were immediately moved to the LPF_TAGRALIA laboratory in Madrid (Spain). Each batch or sample unit was divided into two subsamples (Fig. 1): Subsample 1 (1 kilogram of fresh olives) was sent to the reference laboratory CM Europa S.L. (Jaén, Spain) to undergo reference analysis of the oil content on a fresh-weight basis using Soxhlet (OCFW SOXHLET (%)); Subsample 2, constituted by 8-9 individual olives, was used for NIR analysis.

In LPF_TAGRALIA facilities, the fresh weight (FW(g)) of each fruit from Subsample 2 was measured using a precision balance (ADP 720/L; Adam Equipment Co. Ltd., Kingston, Milton Keynes, UK). Next, the fruits were dried in an oven (Conterm Poupinel; JP SELECTA S.A., Abrera, Barcelona, Spain) at 105 °C until a steady weight was achieved. Dry whole individual olive fruits were measured using an NMR Minispec NMS100 (Bruker Optik GmbH). Measurements were made in 30-mm-diameter glass tubes. The device was calibrated with 9 quantities of olive oil (from 0.05 g to 1.35 g), and a calibration line was built between NMR responses and corresponding oil weights. The oil content was given directly by software for each fruit in grams (OC NMR (g)). OCFW NMR (%) determined as a percentage on a fresh-weight basis by NMR for every fruit (j) belonging to the batch (i) was computed considering the FW(g) of each fruit according to (1).

= ⋅ OCFW OC g FW g (%) ( ) ( ) 100 (1) i j NMR i,j NMR i,j ,
In total, 278 dried olives from Subsamples 2 were stored in a dark and fresh place until the end of the harvest season. They were then moved to IRSTEA (Montpellier, France) facilities for spectral analysis. The relative reflectance hyperspectral images of each dry olive fruit were acquired using a vision system comprising a line-scan push broom camera (model HySpex SWIR-320m-e; Norsk Elektro Optikk, Skedsmokorse, Norway). The spectral range of the camera was 1000-2500 nm with spectral sampling every 6 nm. However, due to the low ratio signal-to-noise in the extreme of the spectra, only the range of 1003-2208 nm was considered (202 spectral bands). A halogen light source was used to illuminate the material to be imaged by the camera. The halogen bulb was switched on 30 min prior to taking a measurement to stabilize the light-source temperature drift and improve the spatial lighting uniformity. Reflectance images were obtained by scaling HSI images using a standard white. Absorbance images (-log 10 ) were computed. The average spectrum of each fruit was considered for further analysis. Savitsky-Golay smoothing and differentiation algorithm were applied to the absorbance spectra: a polynomial of order three was fitted to a width of 21 wavelengths, and the first derivative function was applied to the smoothed absorbance spectra. estimation models. To estimate the OCFW (%) and OC (g), partial least squares (PLS) regression was applied to the corresponding average spectra. The goodness of each estimation model was evaluated through the coefficient of determination (r 2 ), standard error of calibration (SEC), standard error of cross validation (SECV), ratio of the prediction to deviation (RPD), number of latent variables (LV) and slope. All data analyses were performed using MATLAB_R2015a (The MathWorks, Natick, USA) and Statistica 13.3 (TIBCO Software Inc., California, USA) software.

In this work, two strategies were used to build models to estimate OCFW (%) (Fig. 2a,b):

• MODEL 1 to estimate OCFW SOXHLET (%) in individual olives. The model used OCFW SOXHLET (%) determined for each of the 31 batches as input in PLS regression. In this case, to match each average spectrum per olive with one reference data, OCFW SOXHLET (%) was replicated by the number of olives of Subsample 2 of each batch. From estimations made by MODEL 1, it is possible to estimate the OC (g) on individual olives by (2).

= ⋅ .   OC g OCFW FW g ( ) (%) 100 ( ) (2)
i j SOXHLET i j SOXHLET i j , ,

• MODEL 2 to estimate OCFW NMR (%) in a homogeneous batch. From MODEL 1 estimations of OCFW SOXHLET (%), the 278 individual olives were re-clustered to build more homogeneous batches. Olives were sorted from OCFW SOXHLET MODEL 1 (%) of 14% to 32%, using steps of 1%. Sixteen homogeneous groups (all of them with a sufficient number of olives) were built. MODEL 2 used OCFW NMR (%) determined for each one of the 16 batches as input in PLS regression. OCFW NMR (%) in each new batch was calculated considering the OC NMR (g) and FW (g) of the individual olives (j) belonging to the same batch (i) according to (3).

• OCFW OC g FW g n n of olives of the batch (%) ( ) ( ) 100, ( 3
) i NMR j n i j NMR j n i j 1 , 1 , = ∑ ∑ = = =
To match each average spectrum per olive with one reference data per batch, the mean spectrum per batch was computed ((4).

SP sp n (4) i j n i j 1 , = ∑ =
Figure 2 shows that MODEL 1 is based on 278 cases or individual olives belonging to 31 "heterogeneous" batches, while MODEL 2 is based on 16 cases or "homogeneous" batches. The contribution of each olive (j) to the heterogeneity (h) of its batch (i) was computed according to (5), for the 31 "heterogeneous" batches (MODEL 1), and ( 6) for the 16 "homogeneous" batches (MODEL 2), equations adapted from Esbensena, K. H. et al. [START_REF] Esbensena | The variographic experiment[END_REF] .

h OCFW OCFW OCFW FW FW (%) ( %) (%) (5) i j het i j SOXHLET i SOXHLET i SOXHLET i j l , , ,  = - ⋅ = - ⋅ h OCFW OCFW OCFW FW FW (%) ( %) (%) (6) i j hom i j SOXHLET i SOXHLET i SOXHLET i j l , , ,    where  OCFW (%) i j SOXHLET ,
is the OCFW (%) estimated by MODEL 1 for each olive (j) belonging to batch i, OCFW (%)

i SOXHLET

is the OCFW(%) determined by Soxhlet in the reference laboratory for batch i, FW i is the mean fresh weight of the olives belonging to batch i, and

 OCFW (%) i SOXHLET
is the mean OCFW(%) estimated by MODEL 1 for batch i according to (7).

  = ∑ = OCFW OCFW n (%) (%) (7)
i SOXHLET

j n i j SOXHLET 1 ,
On the other hand, two models were computed to estimate OC (g) in individual olives (Fig. 2c,d): • MODEL 3 to estimate OC SOXHLET (g) in individual olives. The model used OC SOXHLET (g) determined for each fruit and calculated according to (8), as input in PLS regression.

= ⋅ OC g OCFW FW g ( ) (%) 100 ( ) (8)
i jSOXHLET iSOXHLET i j , ,

• MODEL 4 to estimate OC NMR (g) in individual olives. The model used OC NMR (g) directly determined by NMR for each fruit as input in PLS regression.

Results and Discussion

Reference analysis. The total range of OC NMR (g) for individual olives varied from a minimum of 0.09 g up to a maximum of 1.29 g (n = 278), that is, all of olives were within the range of the calibration curve (Fig. 3). Such an OCNMR range is even wider than that reported by de la Rosa, et al. [START_REF] De La Rosa | Fruit characteristics and fatty acid composition in advanced olive breeding selections along the ripening period[END_REF] (0.1 g to 0.9 g in individual fruits). Table 1 shows that the average olive FW was 2.07 g (±0.74) similar to the average FW per fruit reported by de la Rosa, et al. [START_REF] De La Rosa | Fruit characteristics and fatty acid composition in advanced olive breeding selections along the ripening period[END_REF] (2.13 g). The degree of heterogeneity of FW, in this olive fruit population, was high, with a CV of 35.88%. The average OC NMR per fruit was 0.51 g (±0.24); de la Rosa, et al. [START_REF] De La Rosa | Fruit characteristics and fatty acid composition in advanced olive breeding selections along the ripening period[END_REF] reported an average OC per fruit of 0.49 g. Also in Table 1, OCFW (%) showed average values 23.77% (±5.23) and 23.96% (±3.48) as determined by NMR and Soxhlet respectively. The CV values being 22% for NMR and 14.54% for Soxhlet. Deblangey, et al. [START_REF] Deblangey | Comparative study of two methods (hexane extraction and NMR) for the determination of oil content in an individual olive fruit[END_REF] report a range of variation for OCFW NMR from 17.1% to 35.5% which is much wider in our study: 10.41-46.59%.

Figure 4(a) shows the correlation between FW (g) and OC NMR (g). A strong uphill linear relationship is found r 2 = 0.87 between these variables; even a stronger relationship (r 2 = 0.94) was found by de la Rosa et al. [START_REF] De La Rosa | Fruit characteristics and fatty acid composition in advanced olive breeding selections along the ripening period[END_REF] in individual olives characterized for a breeding program. As expected, the quantity of oil (absolute value -g-) is higher for larger fruits, while OCFW NMR refers to oil concentration.

Figure 4(b) shows the correlation between FW(g) and OCFW NMR (%) of individual olives. The determination coefficient (r 2 ) is equal to 0.25, indicating that the relationship between the FW(g) of an individual and its OCFW(%) is moderate but non-relevant. The work of de la Rosa, et al. [START_REF] De La Rosa | Fruit characteristics and fatty acid composition in advanced olive breeding selections along the ripening period[END_REF] confirms this result with an r 2 of 0.005 when comparing FW (g) of one fruit vs the OCFW (%) of its batch.

oil content referred to fresh weight in heterogeneous batches. Using spectral data, a first PLS model (MODEL 1) was built for estimating the oil concentration (OCFW (%)) since the OCFW (%) determined by batch is the most usual information used by industry to characterize the product previous to processing. The inputs for MODEL 1 covers the spectral information of each individual fruit plus one common reference per batch, OCFW SOXHLET (%). Figure 5 shows the parameters that characterize the performance of MODEL 1. The coefficient of determination is low with r 2 = 0.67, that is, the explained variance of the model is only 67% even when the number of LV is high (12). SECV is 2% and the RPD is between 1.5 and 2, indicating that the model can only segregate between high and low values [START_REF] Williams | Tutorial: Items to be included in a report on a near infrared spectroscopy project[END_REF] , which agrees with previous works.

Figure 5 shows the actual values of OCFW SOXHLET (%) as compared to those estimated by MODEL 1. A vertical dispersion of data is found for the olives belonging to the same batch. The intra-batch SD estimated for OCFW SOXHLET MODEL 1 was 1.54% while the inter-batch SD reached 7.98%. Therefore, the intra-batch variability is 19.3% of the inter-batch variability, providing an idea of the heterogeneity of the estimated OCFW SOXHLET MODEL 1 per olive within each batch. A first question arises: is this variability due to an estimation error in MODEL 1 or to the intrinsic heterogeneity between the fruits included in the same batch? Assuming that the NMR technique is the best way to determine the OCFW (%) for individual olive fruits, the answer to this question could be found by comparing the estimations of MODEL 1 with the NMR values for individual fruits.

Figure 6(a) shows the OCFW SOXHLET (%) per batch (31 values repeated to match 278) vs OCFW NMR (%) per fruit (n = 278), pointing to a low relationship (r 2 = 0.44) among them. An important difference was noted in the range of OCFW (%) determined by Soxhlet (18.9-31.5%) and that determined by NMR (10-47%), evidenced by the low value of the slope (0.44) of the fitted line (1 for the bisector). The intra-batch SD of OCFW NMR was 3.76%, while the inter-batch SD was 11.6%, that is, the intra-batch variability of the OCFW NMR was 32.4% of that inter-batch. Therefore the intra-batch variability is even higher when OCFW (%) is determined by NMR than when estimated by MODEL 1 using Soxhlet as reference analysis. This result is supported by the work of Deblangey, et al. [START_REF] Deblangey | Comparative study of two methods (hexane extraction and NMR) for the determination of oil content in an individual olive fruit[END_REF] , which, under similar conditions, reported the higher sensitivity of the NMR methodology as compared to other reference analysis when determining OCFW (%) for individual fruits. Furthermore, Deblangey, et al. [START_REF] Deblangey | Comparative study of two methods (hexane extraction and NMR) for the determination of oil content in an individual olive fruit[END_REF] established that NMR generates the lowest precision errors. Figure 6(b) shows OCFW NMR (%) per fruit vs OCFW SOXHLET MODEL 1 (%) estimated per olive according MODEL 1; also in this case the batch effect is strongly attenuated (slope of 0.45). The correlation between the estimated values in MODEL 1 and true value was improved (r 2 = 0.61), indicating that the estimations of MODEL 1 are nearer to the actual olive value of OCFW (%) for each olive as compared to that of the OCFW SOXHLET (%) determined per batch. Estimations with MODEL 1 expand the limits of OCFW to a range from 16% to 32%, though the model cannot properly estimate the OCFW (%) for individual olives beyond these limits, leading to saturated estimations especially at its upper limit (46.59%) according to OCFW NMR . Similar limits of OCFW (%) have been found by other researchers, with minimum and maximum OCFW values between 5% and 44% [START_REF] Leon-Moreno | Usefulness of portable near infrared spectroscopy in olive breeding programs[END_REF][START_REF] Salguero-Chaparro | Near infrared spectroscopy (NIRS) for on-line determination of quality parameters in intact olives[END_REF][START_REF] Cayuela | NIR prediction of fruit moisture, free acidity and oil content in intact olives[END_REF] .

Figure 7(a) identifies nine outliers (red circles) when comparing OCFW SOXHLET MODEL 1 vs OCFW NMR (%). In Fig. 7(b), the values of FW (g) vs OC NMR (g) are plotted highlighting the outliers. The absolute values in grams of FW and OC determined for each of the highlighted data are within the range of calibration of NMR, as well as within the weight range for this population, however when combined in the computation of OCFW NMR (%) lead to abnormal values, either being low or high. Considering that the manipulation of the sample is minimal (fresh whole fruit without pretreatment that is weighed in a balance with a scale accuracy of 0.001 g and then directly measured by NMR), it seems that these highlighted values are singular individuals detected by NMR and not measurements errors. Thus, the high heterogeneity in the OCFW NMR (%) of olives is demonstrated to be due to intrinsic differences among fruits even when they belong to the same batch, with MODEL 1 partially detecting such intrinsic differences.

When relating the OCFW SOXHLET (%) of olives estimated by MODEL 1 with the actual value determined by NMR for each olive excluding the singular individuals a determination coefficient of r 2 = 0.69 is obtained improving the r 2 of 0.61 in Fig. 6(b). This fact corroborates the higher accuracy of OCFW (%) estimates with MODEL 1 for each olive with respect to OCFW SOXHLET .

oil content referred to fresh weight in homogeneous batches. To generate homogeneous batches, the fruits were clustered into groups according to the values of OCFW SOXHLET MODEL 1 (%). Figure 8 plots the contribution of each olive to the heterogeneity (h) of its batch according to Eqs 5 and 6. The blue line points a high contribution of individuals to the heterogeneity of the batches when 31 cluster are considered. In the case of considering 16 groups, the heterogeneity (indicated by the red line) stays around 0, and thus selected as best option.

Figure 9(a) shows the average estimate OCFW SOXHLET (%) with MODEL 1 for each cluster (Eq. 7) compared to that of OCFW NMR (%) (Eq. 3) achieving a high determination coefficient (r 2 = 0.97). This means that 97% of the variance of OCFW SOXHLET MODEL 1 is explained by the actual OCFW (%) per cluster determined by NMR. Figure 9(b) shows the performance of MODEL 2, adjusted on the 16 homogenised spectra using OCFW NMR as dependent variable. In this case the coefficient of determination was r 2 = 0.96, with only 3 LV, indicating the high robustness of the model. SECV was 1.2%, RPD was 4.74, and the slope was 0.92, indicating that quantitative predictions are possible even at the extremes [START_REF] Williams | Tutorial: Items to be included in a report on a near infrared spectroscopy project[END_REF] .

Figure 10(b) shows the average treated spectra of the 16 clusters considered as homogeneous, together with the loading of the wavelengths in MODEL 2 (Fig. 10a). Considering that 1200 nm is the main absorption band for fats and oils, a spectral zoom between the positive peak located at 1153 nm and negative peak located at 1231 nm is presented (Fig. 10c). It can be observed that the average spectra of each cluster are ordered according from lower to higher OCFW (%). This is not the case when considering the original batches (data not shown). calibration models to estimate the oil content (g). Figure 11 shows, from left to right, two models calibrated according to Eq. 8 with using OC SOXHLET (g) (MODEL 3) and to OC NMR (g) (MODEL 4) as dependent variables; both figures show a non-linear behaviour (banana-shaped distribution). High OC (g), usually corresponding to the largest fruits, are not accurately estimated with the PLS models and seem saturated above 1 gram of OC per fruit.

Figure 12 confronts the analysis of the residuals in MODEL 4 by means of comparing the actual OC NMR values (g) and the residuals of estimates. In this Figure the histograms of actual OC NMR (g) and of estimate residuals are combined with corresponding scatterplot. Actual OC NMR (g) does not follow a normal distribution; the distribution is positively skewed (skewness = 0.73) with a high occurrence for low OC (g). Thus, a heteroscedastic error may be inferred. Such lack of compliance with the restrictions for a linear regression can justify the low quality of the estimation with PLS models.

As stated by Beer´s law, which is valid only for transparent homogeneous materials, and the more practical approach of the Kubelka-Munk equation: f(C) = log(1/R), the information present in an NIR spectrum is related to the concentration of a given substance in a sample [START_REF] Pasquini | Near Infrared Spectroscopy: fundamentals, practical aspects and analytical applications[END_REF] . Most of the quantitative applications are targeted to determine major constituents in the sample, with usual detection limits of approximately 0.1% (m/m). The interactions of the light with the sample are limited to a restricted volume, implying that the change in the signal intensity is due to the major constituents that are inside this volume, representing the % in mass and not the total quantity of this constituents in the sample. Because spectroscopy is sensitive to concentration, in this case, to OCFW (%), but not to the total quantity of one compound (OC (g)), it could be considered a methodological error to calibrate a model directly with oil quantity in grams.

However, from a practical point of view, it is interesting to assess the mass of oil that is entering a mill. Therefore, it is necessary to answer whether it is possible to satisfy this requirement using NIR spectra. Comparing the OC SOXHLET MODEL 1 (g) (Eq. 2) with the true value determined per olive by NMR (OC NMR (g)) were found a high determination coefficient (r 2 = 0.93) and a low standard error of prediction (0.061 g). Therefore, when gathering the OCFW SOXHLET MODEL 1 (%) estimation and the FW (g) of each fruit in Eq. 2 the best estimation of the oil content in grams for individual olives is obtained. Currently, vision machines are being developed to classify olives according to different quality parameters, such as colour and defects, previous to milling. These machines use multispectral vision cameras to extract the parameters included in the classification algorithms. In this state of development, it is perfectly possible to use the images to estimate the FW (g) of each olive [START_REF] Williams | Tutorial: Items to be included in a report on a near infrared spectroscopy project[END_REF] . As proven above, gathering the OCFW (%) information of the batch (i.e., by Soxhlet) together with the fresh weight (g) of olives would allow the estimation of the OC (g) per fruit, and thus the mas of oil (kg) which is entering the mill. Moreover, the use of multispectral cameras focused on the appropriate wavelengths, will lead to spectral models for OCFW (%) quantification that can be implemented in real-time, providing even a more accurate estimation.

conclusions

The complete seasonal heterogeneity in the OC of a commercial mill was characterized through systematic and representative sampling according to a factory standard (Soxhlet, %), together with the NMR oil actual value assessed on individual fruits, ranging from 10.4 to 45.6%, 0.09 to 1.29 g per olive.

The OC (g) estimated using a laboratory-top NMR instrument with specific calibration, is taken as the actual value in this study with a precision level of 0.017 g (0.8% of FW for an average fruit).

Current batch based assessment of the OC (Soxhlet, %) in mills only reproduces 44% of the underlying heterogeneity, despite being the factory standard used for payment to the farmer.

A PLS spectrometry model (1003-2208 nm) based on individual olives to estimate the OC reproduces 67% of batch variance and 60% of underlying heterogeneity. Therefore, spectrometry on individual olives helps to assess the variability of the oil content (%) in-mill even using batch values as the dependent variable.

It has been corroborated that it is a methodological error to develop PLS spectrometric models to directly estimate the OC (g) of the fruits since spectroscopy is sensitive to concentration but not to the total quantity of one compound (OC, g). However, the estimation of the OC (% fresh weight) by spectrometry on individual olives together with the assessment of fruit fresh weights (g) reproduces 93% of the variance of the oil content (g) in individual olives. A standard error of prediction of 0.061 g in the OC (g) (2.9% of FW for an average fruit) was reached through the combination of spectrometry and weight in individual olives, a value that approaches the potential limit according to the NMR reference (0.017 g) taken as the actual value.

The improvement in dealing with sample heterogeneity provided by the combination of spectrometry and olive fresh weights contributes to the fair rating of the product value, as well as to provide more accurate process settings in the mills.

It may be foreseen that developing olive grading lines combining spectrometric and physical properties of individual olives will become a commercial target in the near future for the olive oil industry.

  optimal management of oil content variability in olive mill batches by niR spectroscopy e. c. correa 1 , J. M. Roger 2 , L. Lleó 1 , n. Hernández-Sánchez 1 , p. Barreiro 1 & B. Diezma 1 total oil content (oc) is one of the main parameters used to characterize the whole of olives entering a commercial mill, quantified by the total fresh weight of the lot and the oil concentration (%) assessed in a representative sample on olive paste, by means of chemical extraction. nuclear magnetic resonance (nMR) and niR spectroscopy are alternative methods even at individual olives. this work evaluates several strategies to calibrate precise niR models for the estimation of the total oc. to this end, 278 olives were analysed covering whole season variability in terms of olive fresh-weight and the corresponding OC by chemical extraction in 31 batches. The average spectra from hyperspectral NIR images (1003-2208 nm) were computed for each fruit and the actual OC (g) of those olives determined by NMR (0.09 to 1.29 g with a precision of 0.017 g). According to the results, current batch based assessment of the OC (Soxhlet, %) in mills only reproduces 44% of the underlying heterogeneity, despite being the factory standard. The incorporation of individual NIR spectra (278) to the 31 Soxhlet values of the batches allows a 67% explanation of the OC (%) of olives. When estimating OC (g) gathering individual fresh weight and the estimation of oil concentration in olives, a standard error of prediction of 0.061 g is reached (r 2 = 0.93), a precision value that approaches the potential limit according to the NMR reference (0.017 g).

Figure 1 .

 1 Figure 1. Scheme of the sampling and methods implemented. FW: fresh Weight (g); OCFW: oil content fresh weight (%); OC: oil content (g); sp: average spectrum per fruit.

Figure 2 .

 2 Figure 2. Schemes of the input and output of each PLS model. Input in MODEL 1 (a), MODEL 3 (c) and MODEL 4 (d) corresponds to 278 cases or olives belonging to 31 batches. Input for MODEL 2 (Eq. 4) corresponds to 16 cases or homogeneous batches of olives. MODEL 2 (b) is calibrated according to the variable determined by Eq. 3. MODEL 3 is calibrated according to the variable determined by Eq. 8.

Figure 3 .

 3 Figure 3. Calibration curve of the NMR instrument. The red line indicates the linear fit to the true values (blue dots) obtained by oil weights. The black dotted lines indicate the total range of OC measured in this population of olives.

Figure 4 .

 4 Figure 4. Scatterplots of FW (g) versus OC NMR (g) at the left (a) and OCFW NMR (%) determined by Eq. 1 at the right (b). The red line indicates the linear fit, and the dotted red lines indicate the confidence levels at 95%.

Figure 5 .

 5 Figure 5. Scatterplot result of PLS MODEL 1. OCFW (%) measured using SOXHLET methods (x-axis) vs OCFW (%) estimated by MODEL 1 (y-axis). The blue line indicates the linear fit, and the dotted blue lines indicate the confidence levels at 95%.

Figure 6 .

 6 Figure 6. Scatter plots of OCFW NMR (%) (Eq. 1) vs OCFW SOXHLET (%) (a) and OCFW SOXHLET predicted by MODEL 1 (b) for individual olives. The red line indicates the linear fit, and the dotted red lines indicate the confidence levels at 95%.

Figure 7 .

 7 Figure 7. (a) Means with error plot using the integer mode to create grouping intervals, displaying outlier data. (b) Scatterplot of FW (g) vs OC NMR (g). Red points indicate the location of outlier data.

Figure 8 .

 8 Figure 8. Contribution of each individual olive to the heterogeneity (h) of its batch. The blue line is h (Eq. 5) when 31 heterogeneous batches are considered, and the red line is h (Eq. 6) in the case of 16 homogeneous batches.

Figure 9 .

 9 Figure 9. (a) Scatterplot of OCFW NMR (%) per batch (Eq. 3) vs OCFW (%) per batch estimated by MODEL 1 (Eq. 7). (b) Scatterplot result of PLS MODEL 2. OCFW NMR (%) per batch (Eq. 3) vs OCFW (%) estimated by MODEL 2 (y-axis) per batch. The red line indicates the linear fit, and the dotted red lines indicate the confidence levels at 95%. Point labels indicate de batch number.

Figure 10 .

 10 Figure 10. (a) b coefficients of PLS MODEL 2 where vertical red lines indicate peaks at 1153 and 1231 nm. (b) treated average absorbance spectra of 16 homogeneous groups plotted in the complete range and (c) it is shown the detail of these average spectra centred in the most informative wavelengths.

Figure 11 .

 11 Figure 11. (a) Scatterplot result of PLS MODEL 3. OC SOXHLET (g) (Eq. 8) vs OC (g) estimations per fruit (y-axis) (b) Scatterplot result of PLS MODEL 4. OC NMR (g) vs OC (g) estimations per fruit (y-axis). The red line indicates the non-linear fit.

Figure 12 .

 12 Figure 12. Scatterplot with marginal histograms of OC NMR (g) vs residuals of MODEL 4 (g). The red line indicates the non-linear fit.

Table 1 .

 1 Main statistics of the analytical parameters.

		Valid n	Mean	Minimum	Maximum	±Std. Dev. CV (%)
	FW (g)	278	2.07	0.77	4.51	0.74	35.88
	OC NMR (g)	278	0.51	0.09	1.29	0.25	49.60
	OC SOXHLET (g) by Eq. 8	278	0.51	0.16	1.33	0.24	46.61
	OCFW NMR (%) by Eq. 1	278	23.77	10.41	46.59	5.23	22.00
	OCFW SOXHLET (%)	31	23.96	18.92	31.54	3.48	14.54

Scientific RepoRtS | (2019) 9:13974 | https://doi.org/10.1038/s41598-019-50342-6 www.nature.com/scientificreports

Scientific RepoRtS | (2019) 9:13974 | https://doi.org/10.1038/s41598-019-50342-6

This approach proves that OCFW(%) estimates for individual fruits with MODEL 1 are accurate enough for classification purposes and can be used to generate homogeneous groups to reconfigure batches for a reference analysis, or either to select olives for breeding purposes.

Acknowledgements

The funding of this work has been covered by Comunidad de Madrid and European Union through S2013/ABI-2747 (TAVS-CM) project and by Universidad Politécnica de Madrid (Mentoring Program). We are thankful for the technical help by David Perez of the mill LA PONTEZUELA SLU (Spain) and the support of the enterprise and staff of MULTISCAN TECHNOLOGIES SL (Spain). We express our gratitude to Prof. María Gómez del Campo from Universidad Politécnica de Madrid for the use of the NMR oil content measurement equipment.

Author Contributions

Each author has made substantial and significant contributions to this work, has approved the submitted version and has agreed both to be personally accountable for the author's own contributions and to ensure that questions related to the accuracy or integrity of any part of the work, even ones in which the author was not personally involved, are appropriately investigated, resolved, and the resolution documented in the literature. Specifically, the contribution of every author has been as follows: Correa, Lleó, Hernández-Sánchez and Diezma: were responsible of experiment design and carried out the experiments and tests. Correa and Diezma carried out and supervised the data analysis and the proper presentation and interpretation of the results. They wrote the main text of this manuscript and prepared the Figures and Tables. Correa put paper together and submitted the paper to the journal. Lleó, Hernández-Sánchez, Roger and Barreiro made intellectual contributions to the data analysis and contributes to data interpretation. Correa, Diezma, Barreiro and Roger reviewed each paper draft. All authors reviewed the manuscript.

Additional Information

Competing Interests: The authors, members of the LPF_TAGRALIA research team, declare as financial support the research support for data collection: the olives were provided by the mill LA PONTEZUELA SLU (Spain) and the generation of initial lots with respect to the external colour of the olives (green, purple and black) were carried out using equipment owned by MULTISCAN TECHNOLOGIES SL (Spain). These authors declare that they have no conflicts of non-financial interests. The author Jean-Michel Roger declares no conflicts of interest, neither financial nor non-financial.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.