domain: spi phys

text:
Hugo Jonquière

Laurent M Mugnier

Renaud Mercier-Ythier

Study of linear phase shift algorithms and application to deflectometry

For a number of optical methods, including deflectometry, information is encoded in fringe patterns. To extract the phase data from intensity values, a number of phase shift algorithms have been designed. In deflectometry, the nonlinearity of the fringe display implies that the displayed fringe pattern presents a harmonic content even if the input pattern is perfectly sinusoidal. The propagation of these harmonics through phase shift algorithms creates parasitic fringe patterns, reminiscent of the initial fringe pattern on the estimated phase. This phenomenon, known as print-through, has been identified as a serious performance limitation.

In this paper, we revisit Surrel's work on harmonic insensitive phase shift algorithms and demonstrate that the class of Discrete Fourier Transform (DFT) phase shift algorithms he defines is very appropriate for the field of deflectometry. We show how to choose the most suitable one depending on the application by performing a complete modeling of the harmonic print-through phenomenon for these DFT algorithms and studying the error propagation for shot noise and temporal perturbations. In a deflectometry context, we demonstrate by means of simulations that carefully chosen DFT algorithms can simultaneously be robust to print-through and perform better with respect to noise than the state of the art nonlinear phase shift algorithms.

Lastly, by comparing experimental mirror shape measurements of the matrix of the secondary mirror of the European Extremely Large Telescope made on the one hand by DFT deflectometry and on the other hand by phase shift interferometry, we demonstrate that the use of DFT algorithms can substantially improve the high spatial frequency measurement capabilities of a deflectometry setup, without the need for a calibration of the display's nonlinearities.

Introduction

The fabrication of an optical part is an iterative process in which shape metrology acts as a feedback. The amplitude and spatial frequencies of the shape defects vary largely during the fabrication of the optical part. To follow the dynamics of the surface, a wide range of metrology processes have been developed, including profilometry, 3D scanning machines, deflectometry, and interferometry [START_REF] Harding | Handbook of optical dimensional metrology[END_REF].

Deflectometry is a non contact, low cost and full field slope measurement process adapted to specular surface. It has been used in optical fabrication as a high spatial frequency measurement tool [START_REF] Huang | X-ray mirror metrology using SCOTS/deflectometry[END_REF][START_REF] Su | SCOTS: A reverse Hartmann test with high dynamic range for Giant Magellan Telescope primary mirror segments[END_REF]. Unlike interferometry, deflectometry measurement does not require a reference wavefront to measure the shape of an optical part, which makes it suited for aspheric or even freeform mirror metrology. In contrast, in interferometry, the reference wavefront is generated by a costly wavefront corrector.

For a number of optical methods, including deflectometry, information is encoded in fringe patterns. To extract the phase data from intensity values, a number of phase shift algorithms have been designed. They rely on several phase shifted fringe patterns and find uses in, for example, phase measuring profilometry [START_REF] Srinivasan | Automated phase-measuring profilometry of 3-D diffuse objects[END_REF], deflectometry [5], and phase shift interferometry (PSI) [START_REF] Johannes | New compensating four-phase algorithm for phase-shift interferometry[END_REF].

In linear phase shift algorithms, phase is calculated as the argument of a complex linear combination of shitfed fringes. Therefore these algorithms allow for a choice of parameters: number of steps, phase steps, coefficients of the linear complex combination. These parameters can be chosen to minimize bias and noise in the estimated phase. The 4-image algorithm has been used in deflectometry setups [START_REF] Zhou | Improved method for rapid shape recovery of large specular surfaces based on phase measuring deflectometry[END_REF][START_REF] Chunfeng | Improved phase-measuring deflectometry for aspheric surfaces test[END_REF][START_REF] Su | Software configurable optical test system: a computerized reverse hartmann test[END_REF], while phase shift interferometry usually uses more images in order to reduce the impact of noise, phase step errors and vibrations [START_REF] Huntley | Suppression of phase errors from vibration in phase-shifting interferometry[END_REF][START_REF] Hibino | Susceptibility of systematic error-compensating algorithms to random noise in phase-shifting interferometry[END_REF][START_REF] De | Derivation of algorithms for phase-shifting interferometry using the concept of a data-sampling window[END_REF].

In deflectometry, as in fringe projection profilometry, the nonlinearity of the fringe display implies that the displayed fringe pattern presents a harmonic content even if the input pattern is perfectly sinusoidal [START_REF] Hack | Measurement uncertainty of linear phase-stepping algorithms[END_REF][START_REF] Servin | Fringe Pattern Analysis for Optical Metrology[END_REF]. The propagation of these harmonics through phase shift algorithms has been identified or defocusing propose to reduce the amplitude of the harmonics before acquisition [START_REF] Peng | Phase error correction for fringe projection profilometry by using constrained cubic spline[END_REF][START_REF] Zeng | Complex surface three-dimensional shape measurement method based on defocused gray code plus phase-shifting[END_REF], but these methods require an extra step in the measurement procedure. To avoid propagation of this print-through while requiring no calibration procedure or defocusing, state of the art algorithm use nonlinear phase shift algorithm with exponential fringes [START_REF] Liu | Gamma model and its analysis for phase measuring profilometry[END_REF][START_REF] Niu | Compensation of nonlinear errors of phase-shifted fringes in phase measuring deflectometry[END_REF][START_REF] Chen | Generic exponential fringe model for alleviating phase error in phase measuring profilometry[END_REF]. This method performs almost as well as the gamma calibrated method with respect to noise and harmonic propagation [START_REF] Babaei | Exponential fringe pattern projection approach to gamma-independent phase computation without calibration for gamma nonlinearity in 3D optical metrology[END_REF].

Another method to limit print-through stands in linear phase shift algorithms designed to be insensitive to a number of harmonics. In particular in this paper, we revisit Surrel's work [START_REF] Surrel | Design of algorithms for phase measurements by the use of phase stepping[END_REF] and demonstrate that the class of Discrete Fourier Transform (DFT) phase shift algorithms he defines is relevant in the field of deflectometry. Using a calculation of the sensitivity of linear Phase Shift algorithms, Surrel details how to design DFT algorithms insensitive to harmonics up to an arbitrary order. To study the specific case of DFT algorithms and demonstrate their relevancy, we offer in this paper a complete modeling of the harmonic print-through phenomenon for DFT algorithms, which includes error propagation for shot noise, temporal perturbations and harmonics.

We then perform with numerical simulations a performance comparison of N-DFT algorithms versus nonlinear phase shift algorithms with respect to noise. Results demonstrate that our above analysis allows us to choose within the family proposed by Surrel the optimal algorithm for a typical deflectometry setup, whose performance is limited by the nonlinear display [START_REF] Jiang | Nonlinear phase error compensation for fringe deflectometry measuring system[END_REF]. We find these algorithms to be superior in terms of phase error residuals to all the state of the art uncalibrated algorithms.

Lastly we apply N-DFT algorithms to experimental deflectometry data and compare them to phase shift interferometry measurements. We demonstrate that choosing the relevant DFT algorithm within the family proposed by Surrel significantly improves the high frequency measurement capacity of our deflectometry setup.

2 Phase Shift algorithm: print-through

Phase Shift algorithm formalism

We define a intensity encoded phase pattern I as a 2π periodic function of a phase φ:

I : R → R φ → I(φ) .
Phase shift algorithms estimate the phase φ from a series of phase shifted intensity patterns I(φ + δ n ) where δ n is the n-th phase shift. These phase steps can be performed through time, but single-shot phase shift methods have also been developed which rely on spatial phase steps [25,[START_REF] Nguyen | Single-shot deflectometry for dynamic 3D surface profile measurement by modified spatial-carrier frequency phase-shifting method[END_REF]. Because it is 2π-periodic, the intensity pattern can be decomposed in a φ dependent Fourier series:

I(φ) = +∞ k=-∞ α k e ikφ . (1)
Linear Phase Shift algorithms aim to estimate φ as the argument of a complex linear combination of N phase shifted intensity patterns:

φ = arg( N -1 n=0 c n I(φ + δ n )). (2)
Using the Fourier decomposition of

I(φ + δ n ): φ = arg( N -1 n=0 ∞ k=-∞ c n α k e ik(φ+δn)) , φ = arg( ∞ k=-∞ α k e ikφ [ N -1 n=0
c n e ikδn ]).

Following Surrel [START_REF] Surrel | Design of algorithms for phase measurements by the use of phase stepping[END_REF], we define the characteristic polynomial:

P = N -1 n=0 c n X n . ( 3
) φ = arg( ∞ k=-∞ α k e ikφ [ N -1 n=0 c n (e ikδ ) n ]), φ = arg( +∞ k=-∞ α k e ikφ [P (e ikδ ]). (4)
This characteristic polynomial of degree N behaves as a frequency transfer function in the spectral plane. Indeed, each initial Fourier coefficient α k is filtered (multiplied) by P (e ikδ ), which does not depend on φ. This characteristic polynomial is fully defined by our choice of complex coefficients c n , n ∈ [1 : N ]. Additionally, if the phase step is a multiple of 2π N , the characteristic polynomial is N-periodic in the frequency plane:

P (e i(k+N )δ ) = N -1 n=0 c n e i2π n(N +k) N = N -1 n=0 c n e i2π nk N = P (e i(kδ+N ) ).
For the estimator φ to retrieve the correct phase, the necessary and sufficient condition is:

∀φ, e iφ ∝ +∞ k=-∞ α k e ikφ [P (e ikδ ].
The Fourier vectors being a basis, the necessary and sufficient condition can be written:

P (e ikδ ) ∈ R if k = 1 = 0 if k = 1 .
To ensure that our estimator retrieves the correct phase in noiseless conditions, the characteristic polynomial should have an infinite number of roots. Since this polynomial is of degree (N-1), some harmonics can propagate as errors. In fringe projection profilometry as in deflectometry, the nonlinear display of fringe creates an important harmonic content, and design of algorithms insensitive to harmonics is therefore relevant [5,[START_REF] Wu | Phase error analysis and reduction in phase measuring deflectometry[END_REF][START_REF] Liu | Gamma model and its analysis for phase measuring profilometry[END_REF][START_REF] Niu | Compensation of nonlinear errors of phase-shifted fringes in phase measuring deflectometry[END_REF]. An algorithm will be insensitive to the j-th harmonic if:

P (e ikδ ) = 0 if k ∈ -j, j, k = 1.
Using polynomial roots, Surrel demonstrated that phase shift algorithms defined by N phase steps with:

c n = e -i2πn N δ = 2π N (5)
are insensitive to any harmonic up to the (N-2)-th. Surrel names these algorithms the Discrete Fourier Transform algorithms (DFT). Since these algorithms are fully defined by the number of intensity patterns N, and in order to distinguish each DFT algorithm, we will refer to them as N-DFT algorithms.

The characteristic polynomial acting as a Frequency transfer function, we can visualize the sensitivity of any given algorithm to the propagation of each harmonic k [START_REF] Servin | Fringe Pattern Analysis for Optical Metrology[END_REF]. Figure 1 shows the absolute value of frequency transfer function H N (k) = P (e ikδ ) of N-DFT algorithms for different values of N in the range k ∈ [-15, 15], N ∈ [START_REF] Srinivasan | Automated phase-measuring profilometry of 3-D diffuse objects[END_REF][START_REF] Chunfeng | Improved phase-measuring deflectometry for aspheric surfaces test[END_REF][START_REF] De | Derivation of algorithms for phase-shifting interferometry using the concept of a data-sampling window[END_REF][START_REF] Manh The Nguyen | Non linearity response correction in phase-shifting deflectometry[END_REF]. As explained above and illustrated on Figure 1, a N-DFT algorithm is sensitive to harmonic k if its frequency transfer function evaluated at k (k = 1) is not zero. In the range k ∈ [-15, 15], the 4-DFT algorithm is sensitive to harmonics [-15,-11,-7,-3,5,9,13], while the 8-DFT algorithm is sensitive to harmonics [-15,-7,9]. In the same range, the 12-DFT algorithm is sensitive to harmonics [-11,13] and the 12-DFT algorithm is sensitive to harmonic -15 only. Every N-DFT algorithm is insensitive to the first (N-2) harmonics, as stated by Surrel and al. [START_REF] Surrel | Design of algorithms for phase measurements by the use of phase stepping[END_REF].

Analysis of phase shift estimation error: impact of a nonlinear display

While the harmonic sensitivity can be determined by the characteristic polynomial, the resulting phase error remains to be calculated. Using the previous formalism, we are able to calculate the propagation of any harmonic through the N-DFT algorithm to obtain the corresponding phase error, i.e, the manifestation in the reconstructed phase of the harmonic propagation through the reconstruction (See calculation in Appendix A). Assuming a Fourier distribution dominated by the fundamental frequency, we obtain a measured phase φ:

φ = φ + +∞ k=-∞,k =1 ( α (kN +1) α 1 ) sin(kN φ) + +∞ k=-∞,k =1 ( α (kN +1) α 1 ) cos(kN φ). (6)
Until this point, no assumption was made on the nature of φ: it can be both a temporal or a spatial distribution. For the rest of this paper, we will explicitly apply this calculus for the temporal deflectometry case, in which φ(x, y) is a static spatial distribution and phase steps are performed through time. We underline that all the results in this paper stays valid for any nature of φ and/or any phase step implementations. In this case:

φ(x, y) = φ(x, y) + +∞ k=-∞,k =1 ( α (kN +1) α 1 ) sin(kN φ(x, y)) + +∞ k=-∞,k =1 ( α (kN +1) α 1
) cos(kN φ(x, y)). [START_REF] Zhou | Improved method for rapid shape recovery of large specular surfaces based on phase measuring deflectometry[END_REF] We find that most of the harmonic content is filtered out by the phase shift algorithm, but that some harmonics propagate as print-through of the input phase field on the final phase, as a frequency multiple of N times the fundamental fringes pattern. This print-through occurs always at these specific frequencies, and are characteristic of harmonic propagation in N-DFT algorithms. The amplitude of the print-through does not depend on N, but only on the initial relative amplitude of the harmonic to the fundamental. For example, from Equation ( 6), the 8-algorithm is sensitive to frequencies k such that:

k ≡ 1 (mod 8), k = 1.
In this particular case, frequencies ... -15, -7, 9, 17, ... will propagate as print though on the phase estimation. Since harmonics correspond to a linear combination of a positive and a negative frequency, the 8-algorithm is therefore sensitive to the 7-th and 9-th harmonics. In the general case, a N-DFT algorithm is sensitive to (N+1)-th and (N-1)-th harmonics (modulo N), respectively through the positive and negative frequency. We underline that the frequency of the print-through is dependent on the chosen N-DFT algorithm and on the fundamental frequency, and note that the amplitude of the print-through does not depend on N. For example, the amplitude of the print-through due to the 9-th harmonic would of phase steps.

Analysis of phase shift estimation error: shot noise

Above section has shown that choosing a N-DFT algorithm in the family exhibited by Surrel can limit harmonic print-through propagation. In a deflectometry of profilometry experimental setup, even in the case of a noiseless display, the measured intensity is contaminated by shot noise. In the high flux regime of a typical deflectometry setup, the electronic noises become negligible and shot noise is the dominant noise source, with a standard deviation proportional to the square root of the mean flux (in photo-electrons).

To identify the impact of the N-DFT algorithm choice on noise propagation, we have calculated the variance of phase estimated by the N-DFT algorithm in the Appendix B. Assuming a high signal to noise ratio, sinusoidal fringes and linear display, we obtain a variance of:

V ( φ) = 2 N I 0 , (8)
where I 0 is the mean intensity value of the fringes per recorded intensity pattern, or frame (in photoelectrons). N the number of frames of the N-DFT algorithm. For a given total integration time, there is different ways to distribute the corresponding photo-electron count N I 0 . For example, one can use a 4-DFT algorithm with 10k photo-electrons in each frame or a 5-DFT algorithm with 8k photo-electrons in each frame. This results underlines that if the total number of photo-electron is constant, the variance of the estimated phase is constant: both choices will perform as well in terms of noise propagation. If the shot noise is dominant, as assumed above, then the choice of the N-DFT algorithm doesn't impact the final variance of the estimated phase.

The variance of the estimated phase does not depend on φ. Thus, for a spatially homogeneous intensity distribution I 0 , the distribution of the estimated phase noise is spatially homogeneous. We underline that the phase noise resulting from the propagation of the shot noise through N-DFT algorithms is a white noise since shot noise is white and neighboring pixels are not involved in the phase computation of any given pixel.

Analysis of phase shift estimation error: temporal phase perturbations

Since the noise variance of the estimated phase does not depend of N for a given total photo-electrons number, increasing N appears a viable option to filter out the maximum number of harmonics. Nevertheless, in the case of phase steps performed through time, vibrations, phase step miscalibrations, or temporal variations of intensity of the display create time dependent phase errors which propagate through the phase shift algorithm [START_REF] Huntley | Suppression of phase errors from vibration in phase-shifting interferometry[END_REF]. Using the formalism developed above, we have calculated the phase estimation error induced by any temporal perturbation.

Let a(t) a temporal phase perturbation. We define A(ν) the continuous Fourier decomposition of a(t):

a(t) = +∞ -∞ A(ν)e i2πνt dν.
Assuming the phase perturbation is small with respect to 1, we obtain (from Appendices C and D):

φ = arg ∞ k=-∞ α k e ikφ P (e ikδ ) + i +∞ -∞ A(ν)P (e i(kδ+2πνk∆ )dν , ( 9
)
where ∆ is the time step associated to the phase step δ via the sampling frequency ν 0 :

∆ = δ 2πν 0 .
Let us define H(k) the generalized frequency transfer function as:

H(k) = P (e i2πkδ ) + i +∞ -∞ B ν P (e i2π(kδ+νk∆ )dν. ( 10
) φ(x, y) = arg( k=-∞ α k e x,y [H(k)]). (11)
Under temporal perturbations, the frequency transfer function is altered and no longer necessarily null. For example, let us calculate the impact of temporal perturbations in the case of perfectly sinusoidal fringes. These fringes corresponds to the γ = 1 case, for which intensity with a temporal phase error is modeled as:

I(φ, t) = I 0 (1 + cos(φ + a(t))).
From appendix D, a N-DFT algorithm would estimate the phase as:

φ(x, y) = φ(x, y) + i × e -i2φ(x,y) +∞ -∞ A(ν)P (e i(-4π N +2πν∆) )dν . ( 12
)
Since the generalized frequency transfer function is no longer null for k = -1, the k = -1 frequency propagates as print-through as twice the fundamental fringe frequency. The amplitude of this printthrough depends on a scalar product between the characteristic polynomial and the spectrum of the phase perturbation.

In the general case (any given γ), the initial harmonic content is not restricted to the frequencies k ∈ {-1, 0, 1}. Since H(k) is not necessarily null, the print-through can occur through an N-DFT algorithm at all harmonic frequency.

Model: error source identification from phase error structure

Table 1 summarizes the different types of errors and their propagation through the N-DFT algorithm.

Dominant error source

Types of phase error for N-DFT algorithm Shot noise Homogeneous white noise Harmonics

Print-through at N times the fringe frequency Temporal perturbation: γ = 1

Print-through at twice the fringe frequency Temporal perturbation: γ = 1

Print-through at any harmonic of the fringe frequency Table 1: Summary: error source vs phase estimated residual.

Since every source of error produces a characteristic residual, the user can easily discriminate error sources and identify the limiting factor of the experimental setup (shot noise, harmonics, or temporal perturbations). If the phase estimation residual is dominated by print-through at N times the fundamental fringe frequency, previous results show that using an N-DFT algorithm with an increased N will lower the phase error.

3 Simulation: performance comparison

State of the art: nonlinear phase shift algorithm

The insensitivity of N-DFT algorithm to the first (N-2) harmonics suggests N-DFT algorithm are interesting candidates for phase estimators for deflectometry or fringe projection profilometry applications. If the first (N-2) harmonics are filtered, the above Section demonstrated that the higher harmonics unfiltered by the generalized frequency transfer function propagate as print-through of the initial fringes. To ensure minimal harmonics propagation, the intensity pattern on the detector should be as close as possible to sinusoidal fringes. Displays used for deflectometry or fringe projection profilometry are nonlinear and characterized by a γ parameter describing the power law between the displayed intensity I disp and the input I Input :

I disp ∝ I γ
Input , where γ typically varies between 2.2 and 3.5. The nonlinearity of the display induces a harmonic content on the detector even for a perfectly linear photo-detector and for perfectly sinusoidal fringe pattern input.

I(φ) = I 0 (1 + cos(φ)) γ = I 0 +∞ k=0 ( γ n ) cos(φ) n ,
where I 0 is an intensity constant and:

( γ n ) = (γ)(γ -1)...(γ -n) n! .
Because cos(φ) n contains (among others) the n-th harmonic of φ, sinusoidal fringes displayed on a gamma nonlinear screen induce a harmonic content. For this reason, several methods including look up table [START_REF] Maa | Investigation of phase error correction for digital sinusoidal phase-shifting fringe projection profilometry[END_REF] and gamma calibration [START_REF] Babaei | Exponential fringe pattern projection approach to gamma-independent phase computation without calibration for gamma nonlinearity in 3D optical metrology[END_REF] have also been developed to ensure the captured intensity pattern on the detector is sinusoidal. These methods require the display to be carefully calibrated.

To bypass the need for calibration, nonlinear phase shift algorithms have been developed [START_REF] Liu | Gamma model and its analysis for phase measuring profilometry[END_REF][START_REF] Niu | Compensation of nonlinear errors of phase-shifted fringes in phase measuring deflectometry[END_REF][START_REF] Chen | Generic exponential fringe model for alleviating phase error in phase measuring profilometry[END_REF]. Instead of limiting the harmonic content on the detector, these methods display intensity patterns adapted to nonlinear displays, namely exponential fringes:

I(φ) = I 0 (e 1+cos(φ) ) γ = I 0 (e γ+γ cos(φ) ).
The phase is then estimated via a nonlinear phase shift algorithm. In particular, Niu and al. [START_REF] Niu | Compensation of nonlinear errors of phase-shifted fringes in phase measuring deflectometry[END_REF] uses a 4 images nonlinear algorithm where phase is estimated as :

φ = arg(log( (I 4 -I 1 )(I 4 -I 3 ) (I 1 -I 2 )(I 3 -I 2 ) ) + i log( (I 4 -I 1 )(I 1 -I 2 ) (I 4 -I 3 )(I 3 -I 2 ) )), (13)
where:

I j = I(φ + j π 2 ) = I 0 (e γ+γ cos(φ+j π 2 ) ).
This estimator retrieves the exact phase in noiseless conditions and Babaei and al. [START_REF] Babaei | Exponential fringe pattern projection approach to gamma-independent phase computation without calibration for gamma nonlinearity in 3D optical metrology[END_REF] indicate that it performs almost as well as gamma calibrated methods for fringe projection with respect to display non-linearity and noise, without the need for a calibration.

Simulation: methodology and conditions

In this section, we assess, by means of simulation, the efficacy of N-DFT algorithms in comparison to state of the art nonlinear Phase Shift algorithm (NLPS). Algorithms studied in this simulation have different numbers of frames. To ensure a fair comparison in terms of noise between algorithms, we consider that the total exposure time for each algorithm is the same regardless of the algorithm, therefore the total number of photo-electrons received in a pixel is the same for each algorithm. For example, for a total budget of 10000 photo-electrons per pixel, the 4-DFT algorithm is simulated with 2500 photo-electrons per pixel, while the 5-DFT algorithm is simulated with 2000 photo-electrons per pixel.

Using a power law display model, we have generated phase shifted fringe patterns both for the N-DFT algorithms and for the NLPS algorithm. The static phase field φ(x, y) is defined on a [960, 960] grid, and varying linearly in the horizontal direction from 0 to 11 × 2π. Intensity patterns for the N-DFT algorithms were modeled as:

I j (x, y) = I 0 (1 + cos(φ(x, y) + j 2π N )) γ + n j (x, y), (14)
where n j (x, y) is an outcome of a shot noise for the j-th frame, approximated as a Gaussian noise of variance equal to the noiseless intensity pattern and I 0 is a constant characterizing the average number of photo-electrons per pixel on the detector. Exponential fringes are modeled as:

I j,exp = I 0,exp (e γ+γ cos(φ(x,y)+j π 2 ) ) + n j (x, y), (15)
where n j (x, y) is similarly an outcome of a shot noise for j-th frame, approximated as a Gaussian noise of variance equal to the noiseless intensity pattern, and I 0,exp a constant characterizing the average number of photo-electrons per pixel on the detector. These fringe patterns and their noise outcomes unwrapping method needs to be implemented. Phases were unwrapped using the algorithm given in [START_REF] Herraez | Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path[END_REF], based on sorting by reliability. After unwrapping, theoretical phase maps φ(x, y) were subtracted from estimated phases and the root-mean-square (RMS) of the residuals computed.

We have studied different Signal to Noise Ratio (SNR) conditions in our simulation: the photoelectron number ranges from 100 to 10 8 . Because the NLPS algorithm is a 4 image algorithm, the maximum number of photo-electrons in each frame has been adjusted between methods in order to compare algorithms using the same total number of photo-electrons conditions. The photo-electron count in the figure 2 corresponds to the maximum total number of photo-electron for a pixel, with every frame of the algorithm taken into account. To cover different cases of screen calibration, simulation results are presented for both γ values of 2.2 and 3.2. In a range of photo-electrons from 10 3 to 10 5 , the phase errors of N-DFT algorithm decrease in accordance with the analytical model: the logarithm of the phase estimation error decreases with a slope of - 1 2 with respect to the logarithm of photo electron count. N-DFT algorithms are noise limited for N ≥ 5 for both value of γ. As photo electron count increases, the variance of the estimated phase decreases until print-through dominates the residual phase error. Since the print-through amplitude does not depend on the photo-electron count, the phase error remains constant as the photo-electron count rises. The 4-DFT algorithm's print-through error takes over noise for a total number of photo-electrons superior to 10 4 for γ = 2.2. In the γ = 3.2 case, print-trough dominate the 4-DFT algorithm phase error for all the simulated range of photo-electrons.

Simulation: analysis of the results

N-DFT algorithms are insensitive to every harmonic up to the (N-2)-th. As the number of frames of the N-DFT algorithm increases, the total amplitude of the print-through decreases. For γ = 2.2, the 5-DFT algorithm phase error is dominated by print-through for a total number of photo-electron greater than 10 6 . Since the print-through amplitude depends on the initial harmonic content, we expect the 5-DFT algorithm to be noise limited on a shorter range for γ = 3.2. Indeed, the phase error is dominated by print-through for a total number of photo-electron greater than 10 5 in the γ = 3.2 case. For N ≥ 7, N-DFT algorithms are noise limited for all the simulated range, for both γ = 2.2 and γ = 3.2.

The phase error for the NLPS algorithm decreases steadily as photo-electron count increases for both γ values, with a slope close to -1 2 . A linear regression gives: ln(σ) = 0.3914 -0.516 ln(N ph ) for γ = 3.2, R 2 = 0.9998 ln(σ) = 0.312 -0.495 ln(N ph ) for γ = 2.2, R 2 = 0.9996

, ( 16
)
where σ is the least square residual, R the coefficient of determination and N ph the photo-electron count.

being sensitive to harmonics, we observe that the phase error of the NLPS algorithm slightly depends on gamma through the shot noise propagation. For a total number of photo-electrons of 10 8 , the phase error of the NLPS algorithm is slightly lower for γ = 2.2 (0.19mrad RMS) than for γ = 3.2 (0.22mrad RMS). We now compare the N-algorithm and the NLPS algorithm performances. The NLPS outperforms the 4-DFT algorithm for a photo-electron count greater than 10 5 for γ = 2.2. As the harmonic content increases, we expect the NLPS algorithm to outperform the 4-DFT algorithm for a lower photo-electron count, and indeed, for γ = 3.2, the NLPS algorithm outperforms the 4-DFT algorithm for a photo-electron count greater than 10 3 . The 5-DFT algorithm and the 6-DFT algorithm are quickly outperformed once the print-through dominates the phase error. Because the 8-DFT algorithm and the 12-DFT algorithm are noise limited for a photo-electron count ranging from 100 to 10 8 , they outperforms the NLPS algorithm on all this range for both γ = 2.2 and γ = 3.2. The simulation demonstrates that the NLPS algorithm exhibit a greater phase error than noise limited N-DFT algorithms, for γ = 2.2 and γ = 3.2.

Overall, the simulation demonstrates that noise limited N-DFT algorithms are superior to the NLPS algorithm, and that N-DFT algorithms are noise limited for N ≥ 8 and γ ≤ 3.2 on a range of photoelectron count from 100 to 10 8 .

An additional effect not taken into account in the above simulation is that in a experimental deflectometry setup, the pattern is out of focus2 . The defocus aberration induces a low pass spatial filter of the intensity pattern which improves the performances of the N-DFT algorithm by filtering out harmonics, but deteriorates the performance of the nonlinear algorithm by altering the exponential fringe Fourier decomposition. The nonlinear response of the detector also has a negative effect on the performance of the nonlinear algorithm as it alters the exponential fringe Fourier decomposition. For all these reasons, we find N-DFT algorithm to be superior to the NLPS algorithm with respect to the phase error, for N ≥ 8, for all realistic photo-electron counts and values of γ.

Usage of N-DFT algorithm in deflectometry: application to experimental data

Deflectometry is full field slope measurement adapted to specular surface. In situ capabilities, fast measurements and low cost makes the deflectometry a promising process in the optical fabrication context, in particular in the high spatial frequency range [5, 2, 3]. A deflectometry measurement extracts the shape of a specular surface from the distortion of known displayed patterns. To carry out a deflectometry measurement, a display system and a detector are necessary. Figure 3 illustrates a deflectometry experimental setup. A screen is placed in front of an optical surface. A camera records the reflected displayed pattern. The detector plane is conjugated with the optical surface. The z axis direction is the direction of the optical axis of the surface. The x and y axis are then defined to form a direct system R.

Principles of deflectometry

Let us consider any specular surface defined by height z(x, y) in a reference frame defined in Figure 3.

A three point association is required to measure the local slope of this surface using Snell-Descartes reflection law (the points are defined in Figure 3):

• The source point S = (x s , y s , z s ) emitting the ray, localized on the display;

• The impact point M = (x m , y m , z m ) of the ray on the mirror;

• The impact point C = (x c , y c , z c ) of the reflected ray on the camera.

Once the geometrical position of these three points is known, one can calculate the local slope of the mirror using the Snell-Descartes law of reflection [START_REF] Huang | X-ray mirror metrology using SCOTS/deflectometry[END_REF]:

∇z(x, y) = ∂z ∂x , ∂z ∂y = F (C, M, S), (17)

Oz Ox

Screen Camera Optical surface To identify which source point emitted the ray associated with the pair (C,M), deflectometry measurements rely on intensity encoded phase patterns and a phase shift algorithm. We underline that despite the fact that the point S is defined by 3 coordinates, adding the constraint that S must be on the surface of the screen allows the point S to be characterized by 2 phase values 4 . Let us consider a screen displaying an intensity encoded phase pattern I(φ). The position of each point S = (x s , y s , z s ) of the display can be characterized by two suitably defined phases φ = (φ x , φ y ). We define f the phase encoding function which associates the phases φ = (φ x , φ y ) to each point S = (x s , y s , z s ) of the display with:

f : R 3 → R 2 S = (x s , y s , z s ) → φ = (φ x , φ y ) f (x s , y s , z s ).
The easiest way to measure these two scalar phases is to successively use two sets of intensity encoded phase patterns. In the following, we consider that the two phases are measured one after the other. We underline that f can be chosen arbitrarily, at the condition that f is a one-to-one mapping. A simple method is to define φ x and φ y as linear functions of the position on the surface of the screen. The input fringes are therefore sinusoidal and have a constant fringe density. We have:

φ x = S,ex,s 2πν f φ y = S,ey,s 2πν f , ( 18
)
where ν f is the fringes frequency and (e x,s , e y,s ) the two vectors characterizing the screen plane surface system as defined by the axis (Oxs, Oys) in 4.

S = f -1 (φ).

Finally:

∇z = F (C, M, S) = F C, g(C), f -1 (φ) . ( 19
)
Adding geometrical information (position of the pupil of the camera and position and orientation of the screen), the slopes values can therefore be estimated for every mirror point. Once the gradient of the surface has been measured, an integration algorithm estimates the optical surface height z(x,y) by solving the inverse problem: ∇z(x, y) → z(x, y).

This problem has been studied for wavefront slopes measured by Shack-Hartmann sensors [START_REF] Southwell | Wave-front estimation from wave-front slope measurements[END_REF][START_REF] Rousset | Wave-front sensors[END_REF] and for deflectometry [START_REF] Hongyu Ren | Least-squares method for data reconstruction from gradient data in deflectometry[END_REF]. Integration can rely either on a modal or on a zonal basis [START_REF] Mochi | Modal wavefront reconstruction from its gradient[END_REF][START_REF] Huang | Zonal wavefront reconstruction in quadrilateral geometry for phase measuring deflectometry[END_REF]. In this study, we use the zonal algorithm in the Southwell configuration [START_REF] Southwell | Wave-front estimation from wave-front slope measurements[END_REF].

4.2 Experimental application: measurement of the matrix of the secondary ELT mirror

Methodology

The secondary Extremely Large Telescope mirror (ELT) is a 4 meter diameter class convex mirror manufactured by the European Southern Observatory (ESO) as a part of the ELT project. This mirror is the largest secondary mirror ever employed on a telescope, and the largest convex mirror ever produced. This mirror is currently being polishing at the site of Safran Reosc (St Pierre du Perray, France). Convex mirrors requires a concave wavefront reference to be measured by means of interferometry, for example a Hindle sphere. To measure the shape of the ELT secondary Mirror, a 2 meter highly aspherical concave blank called matrix5 has been manufactured, and is currently polished by Safran Reosc as a reference wavefront. This mirror is composed of Zerodur and was not coated at the time of the measurements.

During the fabrication process of the matrix (high frequency shape error above the specification), three measurement have been carried out: one phase shift interferometry measurement using a computer generated hologram as the aspheric reference wavefront, and two deflectometry measurements, which do not need any wavefront corrector. We underline that at the time of the measurements the matrix was in an early fabrication stage and that the state of the mirror is nowhere near the final shape quality of the mirror.

The above section has demonstrated that the 8-DFT algorithm is an improvement over the 4-DFT algorithm with respect to print-through within any photo-electron count greater than 4.10 4 . To demonstrate experimentally this result, the mirror has been measured by deflectometry with a 4-DFT algorithm and a 8-DFT algorithm. The display used is a 1280x1024, 19 inches Terra MF190D-L03 which has not been calibrated for gamma. The manufacturer indicates a gamma value of 2.2. The camera used is a FLIR: CM3-U3-13S2M-CS: 1.3 MP, 30 FPS, SONY ICX445, MONO. The objective is a 1/1.8 10-40mm varifocal objective from Edmund Optics. The shutter time is set to 500ms. The photo-metric data is a mean of 15 measurements per frame to limit shot noise. The total photo-electron count per pixel is approximately 1.10 5 for the 8-DFT algorithm and 5.10 4 for the 4-DFT algorithm. Phases (φ x , φ y ) have been defined as linear functions of the position on the display in accordance with Equation 16. The fringe frequencies have been set to (ν f,x , ν f,y ) = (28.6, 16.1)m -1 .

Results and analysis

Figure 4 presents the result of (a) the 4-DFT and (b) the 8-DFT phase measurement (φ x , φ y ) (rad) of the matrix of the secondary mirror the ELT. In order to visualize high spatial frequency defects and potential print-through, a phase map filtered out of the first ten Legendre's polynomials (L10) is presented under each phase map in 2b. Fringe patterns can be identified in the (φ x , φ y ) L10 substracted phase maps retrieved by the 4-DFT algorithm. A measure of the fringe pattern frequencies indicates that the pattern frequencies is at 4 times the fundamental fringe frequency. In accordance with the model of the Section 2, we identify these fringes as print-through of the initial fringe patterns. This print-through is absent on the L10 substracted phase maps retrieved by the 8-DFT algorithms. As simulated in the above section, We now compare the performances of N-algorithms and NLPS algorithm with respect to shape estimation. Figure 5 frequencies and print-through. To compare the different shape estimation methods, difference between the 3 shape maps are presented in the second row of the Figure 5. The interferometric shape estimation has been spatially registered with the deflectometry measurements using the SURF algorithm [START_REF] Baya | Speeded-up robust features (SURF)[END_REF]. Registration is based on a similarity transformation. The 8-DFT algorithm leads to a shape reconstruction in good accordance with the phase shift interferometry measurement: the difference has a RMS value of 34 nm for a RMS value of each map higher than 125 nm. Additionally, a careful examination of the difference between the two maps shows that a substantial part of it is due to an imperfect registration of the two maps, and also includes parasitic fringe patterns on the interferometric measurement.

However, the 4-DFT algorithm shape reconstruction is dominated by a grid: 274 nm RMS residuals between the 4-DFT algorithm measurement and the 8-DFT algorithm measurement and 276.7 nm RMS residuals between the 4-DFT algorithm measurement and interferometry measurement. Again, a measure of the grid frequency demonstrates that the grid frequency is at 4 times the fundamental fringe frequency, in accordance with the model of Section 2. Therefore, the grid we observe on the 4-DFT reconstruction is due to the superposition of the φ x and φ y print-through propagating through integration. This grid is absent in the 8-DFT reconstruction, as we can observe in the difference map between the 4-DFT algorithm and the 8-DFT algorithm. This result confirms that by choosing carefully the N-DFT algorithm, high spatial frequencies measurement capacities of a deflectometry setup can be substantially improved.

Conclusion

In this paper, we revisited the work by Surrel [START_REF] Surrel | Design of algorithms for phase measurements by the use of phase stepping[END_REF] on a family of phase shift algorithms robust to harmonics, namely N-DFT algorithms. In particular, we fully modeled the propagation of harmonics due to nonlinearities, shot noise, and temporal perturbations through these algorithms. We demonstrated that temporal perturbations and harmonics propagate through N-DFT algorithms as print-through, a periodical signal whose frequency is a multiple of the initial fringe pattern frequency. This modeling allows users to discriminate error sources and identify the limiting factor(s) of their experimental setup. Additionally, we showed that the print-through frequency can be pushed away to an arbitrary high band of spatial frequencies by a judicious choice of the N-DFT algorithm and of the fringe frequency ν f . In a deflectometry context, we demonstrated by means of simulations that carefully chosen N-DFT algorithms can simultaneously be robust to print-through and perform better with respect to noise than the state of the art nonlinear phase shift algorithm. Lastly, by comparing mirror shape measurements of the matrix of the secondary mirror of the European Extremely Large Telescope made on the one hand by N-DFT deflectometry and on the other hand by phase shift interferometry, we demonstrated that the use of N-DFT algorithms can substantially improve the high spatial frequency measurement capabilities measurements was found to be 34 nm for an initial shape estimated to be greater than 125 nm RMS.

The quality of the deflectometric measurement using tailored N-DFT algorithms suggests the possibility for interferometry-free mirror metrology in the future, as the low spatial frequencies of the mirror could be retrieved by mechanical methods and high spatial frequencies could be retrieved by deflectometry, without the need for the costly and specific reference wavefront required for an interferometric measurement.

The propagation of noise through phase shift algorithms has been calculated in the case of a stationary noise, in particular in [START_REF] Servin | Fringe Pattern Analysis for Optical Metrology[END_REF], where Servin et al. calculate the variance of the estimated phase through the use of the noise's power spectral density (PSD). Yet this calculation is not applicable here, because shot noise is not a stationary noise, and therefore a PSD can not be defined in this case [START_REF] Papoulis | Signal analysis[END_REF]Sect 9.3]. The propagation of shot noise through N-DFT algorithms is calculated below.

Intensity and noise are modeled as:

I(φ) = I 0 (1 + cos(φ)) + n(φ),
where n j (x, y) is an outcome of noise and I 0 is a constant characterizing the average number of photoelectrons per pixel on the detector. This model leads to, using equation 2:

φ = arg( N -1 n=0 c n (I(φ + δ n )) + n(φ + δ n )).
Additionally, from [START_REF] Surrel | Design of algorithms for phase measurements by the use of phase stepping[END_REF], for a N-DFT algorithm c n = e -i2πn N

, which leads to:

N -1 n=0 c n (I(φ + δ n )) = N -1 n=0 e -i2πn N I 0 1 + e i(φ+δn) 2 + e -i(φ+δn) 2 + N -1 n=0 c n n(φ + δ n ).
For N-DFT algorithms, δ n = 2π N , which leads to:

N -1 n=0 c n (I(φ + δ n )) = N -1 n=0 I 0 e -i2πn N + e iφ 2 + e -i(φ+ 4π N ) 2 + N -1 n=0 c n n(φ + 2π N ).
And using the result that the sum of the roots of unity is 0: Hence:

φ = φ + tan -1 N -1 n=0 2 N I0 n(φ + δ n ) sin(-φ -2πn N ) 1 + N -1 n=0 2 N I0 n(φ + δ n ) cos(-φ -2πn N )
.

For a high signal to noise ratio hypothesis, we have:

|n(φ + δ n ))| I(φ + nδ) ≤ 2I 0
Therefore, assuming a high signal to noise ratio, the first term of the Taylor decomposition gives:

φ φ - N -1 n=0 2 N I 0 n(φ + δ n ) sin(φ + 2πn N ).
Assuming that photon noise is the dominant noise source, i.e., that the electronic noises are negligible, the noise variance is given by the mean number of photo-electron per pixel:

V (n(φ)) = I 0 (1 + cos(φ)). Elementary trigonometry calculation simplifies the sum as:

N -1 n=0 1 N cos 2 φ + 2πn N 2 sin 2 φ + 2πn N = 1 4 .
Finally the variance of the estimated phase is:

V ( φ) = 2 N I 0 .
We underline that the variance of the estimated phase does not depend on N under the assumption of a constant photo-electron budget N I 0 , which is the assumption made in this paper.

C Print-through propagation of temporal perturbation in the general case

Any given temporal perturbation on the phase can be described by a continuous Fourier transform:

a(t) = +∞ -∞ A(ν)e i2πνt dν.
Let ∆ the time step associated to the phase step δ and the sampling frequency ν 0 , we obtain:

∆ = δ 2πν 0 .
Then from Equation 4 we obtain: .

φ = arg
Given the phase perturbation is small with respect to 1:

+∞ -∞ | A(ν) | 1,
we obtain at the first order of perturbation: c n e i(knδ+2πνkn∆) dν .

φ = arg
Using the characteristic polynomial: A(ν)P (e i(kδ+2πνk∆ )dν .

φ = arg
For the γ = 1, intensity with a temporal phase error is modeled as: I(φ, t) = I 0 (1 + cos(φ + a(t)) = α k e ik(φ+a(t)) .

Then from Equation 4we obtain:

φ = arg 1 k=-1
α k e ikφ P (e ikδ ) + i +∞ -∞

A(ν)P (e i(kδ+2πνk∆ )dν .

We define the generalized frequency transfer function of the algorithm (which depends on phase error Fourier decomposition) as:

H(k) = P (e ikδ ) + i +∞ -∞
A(ν)P (e i(kδ+2πνk∆ )dν.

In particular: A(ν)P (e -i 4π N +i2πν∆) dν.

H(1) = N e iφ + i
From Appendix A, we apply the formula for the argument of a Fourier serie dominated by the fundamental frequency:

φ = φ + 0 k=-1 ( z k z 1
) sin((k -1)φ) + ( z k z 1 ) cos((k -1)φ).

With z k = α k H(k), we obtain: φ = φ + (H(-1)) sin(-2φ) + (H(-1) cos(-2φ), φ = φ + i × e -i2φ(x,y) +∞ -∞

A(ν)P (e i(-4π N +2πν∆) )dν .

The output of this filter propagates as print-through at 2 times the fundamental fringe frequency.

Figure 1 :

 1 Figure 1: Frequency transfer function of N-DFT (N=4, 8, 12, 16) vs harmonic number.

Figure 2

 2 Figure 2 presents the result of the simulations and shows the phase RMS error as a function of photo electron count for various N-DFT algorithms and for the NLPS algorithm, for γ = 2.2 and γ = 3.2. Simulations are average values of 20 noise realizations.

Figure 2 :

 2 Figure 2: Phase error (rad RMS) vs photo-electron count for a) γ = 2.2 and b) γ = 3.2

Figure 3 : 3 :F

 33 Figure 3: Deflectometry experimental setup and coordinates definition

Figure 4 :

 4 Figure 4: (φ x , φ y ) (in radians) extracted of the measurement of the matrix of the secondary ELT mirror. Phase shift algorithm used is a) 4-DFT b) 8-DFT. Lower line displays (φ x , φ y ) with the first ten Legendre's polynomials subtracted in order to visualize low amplitude defects and print-through.

  presents the result of (a) the 4-DFT algorithm, (b) the 8-DFT algorithm and (c) phase shift interferometry shape measurement of the matrix of the secondary ELT mirror. The 4-DFT algorithm and 8-DFT algorithm measurement correspond to the processed phases presented in Figure 4. The first 36-th Legendre's polynomials (L36) have been substracted on each measurement to visualize high spatial

Figure 5 :

 5 Figure 5: Shape (in nm) measured of the matrix of the secondary ELT mirror by 4-DFT Deflectometry, 8-DFT Deflectometry, and phase shift interferometry. First 36-th Legendre's Polynomials have been substracted.

2 N - 1 n=0 2 N - 1 n=0

 2121 c n n(φ + δ n ) N I 0 e iφ , φ = arg N I 0 2 e iφ + arg 1 + c n n(φ + δ n ) N I 0 e iφ , φ = φ + arg 1 + 2 + δ n )) N I 0 e iφ .

c

  n e i(knδ+ +∞ -∞ A(ν)e i2πνkn∆ dν)

∞

  k=-∞ α k e ikφ P (e ikδ ) + i ∞ k=-∞ α k e ikφ +∞ -∞A(ν)P (e i(kδ+2πνk∆ dν .Finally:φ = arg ∞ k=-∞ α k e ikφ P (e ikδ ) + i +∞ -∞A(ν)P (e i(kδ+2πνk∆ )dν .Starting from the final equation of the Appendix C:φ = arg ∞ k=-∞α k e ikφ P (e ikδ ) + i +∞ -∞

cos(φ) = 1 is out of the radius of convergence

In fringe projection profilometry, the detector is optically conjugated with the display

Using a calculation based on tangent half-angle formula

For example, for a planar screen model, a display point is characterized by two coordinate positions on the surface of the screen.

One of the most aspheric large mirror ever polished.

Acknowledgements

We thank ESO for supporting us during the polishing of the ELT M2 Matrix and Vincent Michau for fruitful discussions.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

A Print-through propagation of harmonic content

Let z a Fourier series of a periodic function:

If z is dominated by the fundamental Fourier component:

Then:

which leads to:

.

The first term of the Taylor decomposition gives:

The characteristic polynomial acts as a frequency transfer function for the initial intensity pattern, and we have: z k = α k P (e ikδ ). And from [START_REF] Surrel | Design of algorithms for phase measurements by the use of phase stepping[END_REF] we have:


Yves Moussallam
email: yves.moussallam@ldeo.columbia.edu

Talfan Barnie

Álvaro Amigo

Karim Kelfoun

Felipe Flores

Luis Franco

Carlos Cardona

Loreto Cordova

Virginia Toloza

Monitoring and forecasting hazards from a slow growing lava dome using aerial imagery, tri-stereo Pleiades-1A/B imagery and PDC numerical simulation

Keywords: structure from motion, volcanic dome, photogrammetry, dome collapse, dome growth, pyroclastic flow

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

unprecedented resolution, allowing for the calculation of discharge rates and effusion rates in near real time. A simple model fit to the distance between the dome and crater rim predicted relatively accurately the arrival of the dome toe at the crater rim and the onset of dome collapse outside the crater. Simulations of the path and extent that potential pyroclastic density currents (PDC) generated by dome collapse would follow showed that PDC were not directly threatening populated areas. Over its life cycle as of August 2019, the dome growth was punctuated by frequent explosions, averaging around 30 per day, one of which generated a minor 600 m long PDC on 13 to 15 of July 2018. There appears to be a positive correlation between explosion frequency and lava dome growth rate suggesting that both explosive and effusive processes can coexist, operating at different timescales but responding to the same driving force. A positive correlation is apparent between dome growth rate and seismic activity such as the frequency of tremor and long-period earthquakes suggesting that these might be used as proxies to estimate effusion rate. Initial lava dome effusion rates of 1730 ± 110 m 3 /day in January 2018 declined to 100 ± 150 m 3 /day in June 2019. These growth rates are extremely slow when compared to other lava domes, about 300 to 600 times slower than the lava domes at Mt Unzen (1992) and Mt. St. Helens (1980).

P a g e | 3

Highlights:

• Ultra-high-resolution DEM and time series of lava dome growth

• Real time forecasting of initiation of dome gravitational instability

• Simulations of PDC from dome collapse show no immediate threat to populated areas • Simultaneity of explosive and effusive activities, positively correlated

• Strong positive correlation between effusion rate and seismic event frequencies

INTRODUCTION

Lava domes form by the slow effusion of typically felsic and viscous magma at the surface and are commonly associated with hazardous phenomena. The most common hazard is the formation of pyroclastic density currents generated by partial or complete collapse of the dome, also termed Merapi-type pyroclastic flows. Examples of these types of flow include the 1991 eruption period at Unzen volcano (e.g., [START_REF] Sato | Crumbling of dacite dome lava and generation of pyroclastic flows at Unzen volcano[END_REF], the 1994 eruption period at

Merapi volcano (e.g., [START_REF] Abdurachman | Nuées ardentes of 22 November 1994 at Merapi volcano, Java, Indonesia[END_REF], the 1996-97 eruption period at Soufriere Hills volcano (e.g., [START_REF] Cole | Pyroclastic flows generated by gravitational instability of the 1996-97 Lava Dome of Soufriere Hills Volcano[END_REF] and the 1998-99 eruption period at Colima volcano (e.g., [START_REF] Saucedo | Emplacement of pyroclastic flows during the 1998-1999 eruption of Volcán de Colima, México[END_REF]. Lava domes can also be associated with explosive eruptions and blast generating Peléan-type pyroclastic flows such as in the 1956 eruption at Bezymianny volcano (e.g., [START_REF] Gorshkov | Gigantic eruption of the volcano bezymianny[END_REF], the 1980 eruption at Mt St. Helens volcano (e.g., [START_REF] Lipman | The 1980 eruptions of Mount St[END_REF], the 1991 eruption period at Unzen volcano (e.g., [START_REF] Sato | Crumbling of dacite dome lava and generation of pyroclastic flows at Unzen volcano[END_REF], the 1997 eruption at Soufriere Hills volcano (e.g., [START_REF] Woods | The explosive decompression of a pressurized volcanic dome: the 26 December 1997 collapse and explosion of Soufrière Hills Volcano, Montserrat[END_REF] and the 2010 eruption at Merapi volcano (e.g., [START_REF] Komorowski | Paroxysmal dome explosion during the Merapi 2010 eruption: Processes and facies relationships of associated high-energy pyroclastic density currents[END_REF]. Monitoring the growth and behaviour of active lava domes is therefore of primal importance to forecast potential explosive events and mitigate their hazards.

In late December 2017 (exact date unknown), after two years of increased activity with frequent eruptions propelled by the emplacement of a small intrusion of magma at shallow level [START_REF] Moussallam | Unrest at the Nevados de Chillán volcanic complex: a failed or yet to unfold magmatic eruption[END_REF], lava reached the surface at the Nevados de Chillán volcanic complex forming a small lava dome (Fig. 1). Subsequent growth of the lava dome was monitored by bimonthly helicopter overflights operated by the Observatorio Volcanológico de Los Andes del Sur (OVDAS), part of the Servicio Nacional de Geología y Minería (SERNAGEOMIN). As helicopter flights became unavailable to continue the monitoring, we tasked the monthly acquisition of Pleiades-1 tri-stereo optical imagery. We used aerial P a g e | 5 photographs taken during observation flights and tri-stereo optical imagery to construct ultrahigh-resolution digital elevation models (DEM) of the summit area and precisely measure the dome growth rate and extruded volume over a period of nineteen months. During this period, we used these measurements to forecast the onset and location of rock fall events. Finally, we performed numerical modelling to predict the extent of PDC that would be generated by a potential dome collapse or large explosion in order to assess hazard to nearby population.

METHODS

In this paper we (i) apply Structure from Motion (SfM) to make measurements of dome volume and geometry from stereo images, (ii) fit Bayesian models to these measurements using Hamiltonian Monte Carlo techniques to estimate effusion rates and predict the start of dome collapse, and (iii) use 2 layer numerical modelling of pyroclastic flows to predict the potential extent of dome collapse. We discuss each of these three techniques in the following sections.

Digital Elevation Model generation

Aerial photogrammetry

Helicopter flights were performed on 21 December 2017, 09 January 2018, 12 January 2018, 16 January 2018, 23 January 2018, 07 February 2018, 22 February 2018,11 March 2018, 4 th April 2018, 18 th April 2018, 16 July 2018 and 9 th August 2019. Flight height was typically around 3400 to 3700 m altitude and circular flight paths were performed around the summit area in order for aerial photograph to capture the lava dome from all angles (Fig. S1).

Cameras were operated manually. Cameras used to take the photograph were equipped with GNSS (GPS of accuracy strictly better than 10m) to provide an approximate position of the camera. 3D models were constructed after each helicopter overflight using the Agisoft Metashape (Agisoft LLC., St. Petersburg, Russia) modelling software (Fig. S1). The Structure from Motion (SfM) software identifies and matches scale invariant features in each images, performs bundle adjustments to refine the camera positions and construct a georeferenced dense point cloud. The georeferenced dense point cloud is then (1) processed with soft-copy triangulation to reconstruct the scene geometry and create a solid 3D mesh (e.g., [START_REF] Burns | Integrating structure-frommotion photogrammetry with geospatial software as a novel technique for quantifying 3D ecological characteristics of coral reefs[END_REF] and (2) used to construct a DEM. On several occasions however, the internal camera GPS did not work and georeferencing was achieved using easily identifiable features in previously reconstructed DEM as Ground Control Points (GCPs). Pixel resolution on the DEMs varied from 10 cm to 1 m (see supplementary table S1 for details).

Tri-stereo satellite optical imagery processing

The PLEIADES constellation consists of two satellites in Low Earth Orbit imaging the Earth at 0.7 m resolution across a 20 km wide swath in the panchromatic band. The satellites are particularly agile, allowing them to acquire multiple images of a target from different angles during a single overpass [START_REF] Gleyzes | Pleiades System Architecture and Main Performances[END_REF]. This stereoscopic facility has been extensively used to measure volcanic deposit volumes (e.g. [START_REF] Bagnardi | High-resolution digital elevation model from tri-stereo Pleiades-1 satellite imagery for lava flow volume estimates at Fogo Volcano[END_REF][START_REF] Di Traglia | The 2014 Effusive Eruption at Stromboli: New Insights from In Situ and Remote-Sensing Measurements[END_REF][START_REF] Carrara | Post-emplacement dynamics of andesitic lava flows at Volcán de Colima, Mexico, revealed by radar and optical remote sensing data[END_REF]Ganci et al., 2019aGanci et al., , 2019b)). Pleiades images were acquired in tristero mode (three views from one overpass) at roughly monthly intervals from October 2018 to June 2019. Images acquired in February and April 2019 had views of the dome obscured by degassing activity and could not be used to generate DEMs. We used the MicMac Structure from Motion software suite [START_REF] Rupnik | MicMac -a free, open-source solution for photogrammetry[END_REF] to create the DEMs, following the Pleiades examples in [START_REF] Pierrot-Deseilligny | Micmac, apero, pastis and other beverages in a nutshell[END_REF], which we summarise briefly here. First the Rational Polynomial Coefficients (RPCs) which are supplied by the data provider and define the polynomial functions that map 3D geographic space to 2D image space were converted to MicMac format using Convert2GenBundle. Then matching P a g e | 9 points are automatically found across all possible pairs of images. This is accomplished using the Tapioca tool in 'All' mode, which uses the Sift+ (Scale Invariant Feature Transform) algorithm [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF] to find tie points. For this step the images were scaled to a width of 10000 (a scaling factor of ~0.5) to save memory. These tie points overdetermine the location and orientation of the imaging systems and so can be used to improve the estimated RPCs using a least squares Levenberg-Marquardt mimization, implemented in the Campari module, improving their accuracy -this is also known as Bundle Block Adjustment. Finally, the DEM is generated from the images using the adjusted RPCs with Ann matching by the Malt tool using the dense matching module UrbanNME.

Vertical precision for both Aerial and Pleiades DEMs were estimated by randomly sampling one hundred points over a stable area (Nuevo crater) and interpolating values at these points for all DEMs, the median height for each sampling point was calculated and subtracted from each height to give a collection of deviations from the median for each point (Fig. S2). The median of all the aerial and Pleiades absolute deviations were then taken, to give the Median Absolute Deviation, or MAD statistic, which came to 1.00 m for Pleiades and 0.476 m for aerial DEMs.

DEM alignment

The aerial and Pleiades DEMs were found to be misaligned with standard topographic models such as Shuttle Radar Topography Mission (SRTM, [START_REF] Farr | The Shuttle Radar Topography Mission[END_REF], as well as misaligned with each other. We aligned our DEMs using the pc_align tool of the NASA Ames Stereo Pipeline software suite [START_REF] Shean | An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery[END_REF], using the Similarity-Point-to-Point Iterative Closest Point (ICP) alignment method. This solves for rotation, translation and scaling to align one DEM with another. The area to be used for alignment was selected by masking out areas where the DEMs showed evidence of change, was noisy, interpolated or otherwise noticeably inaccurate (Fig. S3). For the dome, we then aligned the 21st December 2017 DEM with SRTM, and aligned all subsequent DEMs with the aligned 21st December DEM. This ensured all our DEMs were aligned with each other as well as being at roughly the right geographic location and scale. In early 2019 it became apparent that a gulley was developing below the dome with a depositional fan at the distal end, which could potentially account for some of the erupted volume. We found that DEMs aligned using the area around the dome were poorly aligned around the gulley and fan, and vice versa, so the DEMs were aligned again using an area surrounding the gulley and fan (Fig. S3). Only data from 16th January 2019 onwards reliably and systematically covered the whole feature, so the 16th January 2019 DEM was chosen as our reference and all subsequent DEMs were aligned with that.

DEM differencing

For the dome, we took the 21st Dec 2017 DEM as our reference surface, above which all volumes are measured. This DEM was then resampled to and subtracted from every other aligned DEM, using the ASP geodiff tool to give a time series of differential DEMs (dDEMs) each with the resolution of the newer DEM. These are effectively thickness maps of the dome

Volume calculations

Dome volumes were calculated by digitising each dome and integrating the dDEMs over these shapes. We estimated uncertainties in these volumes by moving the plan shapes of each dome to a neighbouring crater that is well resolved in all the DEMs (Arrau crater, in Figure 1b), but doesn't show any appreciable topographic change, and measuring the volume for all dDEMs. This gives a number of estimates of a 'dome' of zero volume for each dome shape, from which we took the standard deviation as the error. For the gulley and fan a similar procedure was performed where we moved the shape delineating the feature to nearby regions where no significant volume change is expected. We measured the net volume of the gulley and fan to try to estimate the volume of fresh volcanic material present. The hope being that the negative contribution of material removed from the gulley is compensated by the net positive contribution of deposition in the fan. In reality, remobilised material is unlikely to pack back down to the same volume, so it is likely that the values would be an overestimatewe discuss this further in the results section.

Time Averaged Discharge Rates

Time Averaged Discharge Rates (TADRs) (defined after [START_REF] Harris | Lava effusion rate definition and measurement: a review[END_REF] were calculated by differencing successive volumes and dividing by the time interval between them.

Uncertainties in TADR were found by simple propagation of errors.

Effusion rates

Effusion rates (or the instantaneous discharge rate) were calculated by fitting a model to the time series of volumes and finding its derivative. In the absence of an obvious physical model to fit to the data we chose a Gaussian Process (or GP), which can be considered a distribution over functions [START_REF] Görtler | A Visual Exploration of Gaussian Processes[END_REF]; an infinite dimensional multivariate normal distribution where each dimension maps to a point in the domain on which we wish to define a function.

We adopt a Bayesian approach to fitting the GP, where we seek to calculate the posterior distribution, ∨ , over the parameters of our model in light of our data, :

P a g e | 12

| = ∨

Where ∨ is the likelihood, is the prior distribution over the parameter space and is effectively a normalising constant that can be ignored. The GP is defined by the kernel function that is used to calculate the covariance matrix of any multinormally distributed subset of variables in the GP, and here we choose the simple squared exponential kernel function:

= -1 2 -
Where is the element of the covariance matrix that specifies the covariance between points and , and and are 'hyperparameters' that specify the variability and length scale of the functions drawn from the GP. Gaussian processes have the property that their derivatives are also GPs and so can be included in the covariance matrix (e.g., Riihimäki

and Vehtari, 2010) provided the kernel function is modified appropriately for the elements where one or both points are gradients:

= -1 2 - -1 - = -1 2 - 1 1 - 1 -
Where is the element of the covariance matrix where the element is a gradient and is that where both the and elements are gradients. We thus have a model parameter space with dimensions of , , , !, !′. We use a non centred model (Betancourt and P a g e | 13 Girolami, 2013;[START_REF] Papaspiliopoulos | A General Framework for the Parametrization of Hierarchical Models[END_REF] as this is computationally easier and can avoid some problems during Monte Carlo sampling. Briefly, in a non-centered GP model we place an independent and unit-normally distributed prior on our volumes ! and gradients ! at our points and transform them using the Cholesky factor of the covariance matrix. The model was fit to standardised data, and the following priors were selected:

#$%&'( 0,1 #+ ,-..-5,5 0 #$%&'( 0,1 = 1 , 2 = 3ℎ$( 516 7 ! !′ 8 = 2. 0
For the likelihood, we assume our measured volumes are distributed about those predicted by the GP by their observed uncertainties as calculated above in section 2.4: ! ,:;< #$%&'( ! , = ,:;< This model was implemented in the Stan probabilistic programming language [START_REF] Carpenter | Stan: A Probabilistic Programming Language[END_REF], which uses Hamiltonian Monte Carlo to draw samples. These samples were then used to calculate expectations over the posterior distribution (or PD) such as the mean and standard deviation at our points of interest , . Additionally, for every sample ! we generate an additional sample ! + #$%&'( 0, = which gives us a draw from the Posterior Predictive Distribution (PPD). One can think of the PD as giving the best estimate of the volumes and effusion rates in light of our priors and data, while PDD gives the measurements we might expect to make given the above, and, additionally, the measurement error.

Predicting the onset of dome collapse

During the first months of extrusive activity, the lava dome was filling and remained contained in a crater formed by the last two years of explosive activity [START_REF] Moussallam | Unrest at the Nevados de Chillán volcanic complex: a failed or yet to unfold magmatic eruption[END_REF]. As long as the dome was contained within this crater there were no external reason for it to become unstable and collapse (i.e. the dome might still become unstable due to internal overpressure but not due to gravity). In April 2018, DEMs produced by aerial photogrammetry and resulting growth rates were used to forecast future dome growth and the time at which the dome was expected to reach the crater rim, hence becoming partially unstable and starting to generate rock falls (Fig. 5). To predict the time at which this was to happen, we developed a simple conceptual model of the dome growth in a geometry constrained by an inverted conical-shaped crater of height H, radius R, volume V and half angle θ (Fig. 5 C). Within the crater, the level of the foot of the dome forms a horizontal surface at height h, of radius r, that forms the top of a smaller inverted cone of volume v. The horizontal distance d between the dome edge and crater rim can be expressed as a function of effusion rate E and time t as:

= ? -@ 3B '# t D E
Using this relationship, we then fitted the data presented in Fig. 5A using a function of the form:

= F -GH -I E
As before the model was fitted using Stan, as above, although in this instance flat priors were used. During the extrusive activity, a possible scenario was that as the lava dome growth continued, larger portions of the dome might eventually become gravitationally unstable and partial or complete dome collapse might occur, associated with pyroclastic density currents (PDC). In order to determine the spatial area likely to be affected by potential PDC we used the twofluids version of the VolcFlow numerical model [START_REF] Kelfoun | A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents[END_REF][START_REF] Kelfoun | Simulation of block-and-ash flows and ash-cloud surges of the 2010 eruption of Merapi volcano with a two-layer model[END_REF] that simulates concentrated part of PDC (i.e. the block-and-ash flow), the dilute part (ash-cloud surge) and their interaction. Surge can be formed from the concentrated part and, inversely, can form concentrated parts by sedimentation. For details on the method and on equation used, the reader is referred to [START_REF] Kelfoun | A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents[END_REF].

Computer simulation of dome-collapse and explosion induced PDC

We simulated two scenarios; a dome collapse and an explosion. For both scenarios, the total volume of rocks forming the concentrated flow was of 6.4×10 6 m 3 which is the size of the Merapi 2010 Eruption. This size was used as to represent a hypothetical worst-case scenario and is not related to the current dome volume which is an order of magnitude smaller.

Although we stress that larger eruptions are of course possible. Concentrated pyroclastic currents were simulated starting from a circular zone of 50 m in radius for the dome collapse scenario, and a zone of 500 m radius for the explosion scenario. These starting locations were used to simulate eruption initiated by dome collapse, hence starting at the current dome location in the first case and eruption from explosion that could happen anywhere within the general summit location in the second case. At the source area, the initial velocity is zero. The concentrated mass then accelerates along the topography forming a surge. Since model results are strongly dependent on the mass rate at the source, two initial durations of genesis were used for each scenario: 60 and 600 seconds. Short durations (i.e., high mass rates) simulations form widespread concentrated flows and powerful surges. Long durations (i.e., low mass rates) simulations form very small surges and flows, channelized in a few drainages valleys P a g e | 16

and reaching long distances. Fourteen parameters must be defined to run the two-fluids simulation. Some parameters are known (e.g. gravity or topography) or can be estimated from field observations, at least within a given range (e.g. density of the PDCs or volume of the dome that collapse). Five parameters are unknown: three rheological parameters that rule the PDCs dynamics, the exchange coefficient from the concentrated to the dilute PDC, and the density of the mixture formed by the fine particles of the concentrated PDC with volcanic or atmospheric gases. These parameters have been estimated by reproducing real emplacements of the 2010 PDCs at Merapi volcano and the 1997 PDCs at Montserrat, and are similar for both fields cases [START_REF] Kelfoun | Simulation of block-and-ash flows and ash-cloud surges of the 2010 eruption of Merapi volcano with a two-layer model[END_REF][START_REF] Gueugneau | Investigation of surge-derived pyroclastic flow formation by numerical modelling of the 25 June 1997 dome collapse at Soufrière Hills Volcano[END_REF]. For the present study, we have chosen to use the parameters of the 2010 eruptions of Merapi volcano.

RESULTS

The time series of aerial and Pleiades DEMs show the growth of the dome from a small mound at the base of the summit crater in December 2017 to overflowing sometime between April and July 2018 (Fig. 2). Cross sections across the DEMs show the dome maintained a relatively flat surface, consistent with an extruding fluid with a finite yield strength (Fig. 5B), until explosions from January 2019 onwards began to excavate a crater on the top of the dome. Once the dome reached the rim of its host crater, a combination of rock falls and pyroclastic flows appears to have cut a gulley into the volcano's northern flank, which is terminated by a depositional fan. These features are visible as negative and positive height changes in the dDEMs in Fig. 3., respectively. Cross sections of the gulley and fan in the DEMs reveal the progressive incision of the former and growth of the latter from March to June 2019. However, it should be noted that snow is present on the lower flanks of the volcano, so changes in elevation may also partially be a function of changes in snow thickness.

P a g e | 20

Dome volumes and their uncertainties are given in Table 1 and shown in Fig. 4. We were able to estimate dome volumes with uncertainties of approximately 4% and 3% for Aerial and Pleiades DEMs, respectively. The dome grew steadily, reaching a maximum of 3.97e5m 3 in August of 2019. A small decrease between May and June 2019 is presumably due to increasing collapses and explosions. These volumes gave Time Averaged Discharge Rates (TADR) values that fluctuate quite substantially, varying between 1.17e3 and 2.39 e3 m 3 day -1 in the first 3 months, before declining over the rest of the observation period. Uncertainties in TADR averaged 33% for Aerial dDEMs and 121% for Pleiades dDEMs, although the high latter value is largely a function of low effusion rates during the Pleiades acquisition period.

By fitting the GP model we were able to estimate the instantaneous effusion rate, which was found to be 1.73e3 m 3 day -1 on the 9 th January of 2018 before declining to around ~ 1e3 m 3 day -1 and remaining roughly constant from early May to July 2018 before declining again more rapidly in mid-2019. Effusion rate estimates using the GP model had uncertainties of between 3% and 422%.

Seismic activity also shows a remarkably similar pattern to the dome effusion rate during the observation period (Fig. 4D). The total number of seismic event decreases sharply during the first three to four months before reaching a plateau with fairly constant activity from April Fig. 5. The positive correlation between the modelled effusion rate and the frequency of longperiod earthquakes, explosions and tremors is apparent (with R 2 of 0.9, 0.7 and 0.5 and Pvalues of 10 -10 , 10 -6 and 10 -4 respectively). Looking in more details at Fig. 5, three different regimes seem to emerge. The first is for effusion rates below ~100 m 3 /day, the number of seismic events (long-period earthquakes, explosions, and tremors) is high in comparison with the effusion rate and increases quickly with increasing effusion rate. The second regime, for effusion rates between ~100 and ~900 m 3 /day shows no correlation between effusion rate and seismic activity with the number of long-period earthquakes, explosions and tremors remaining fairly constant over this large range of effusion rates. The third regime, for effusion rates between ~900 and ~1700 m 3 /day shows again a strong positive correlation between effusion rates and numbers of explosions, tremors, and long-period earthquakes but with an increase in seismic event with increasing effusion rate that is slower from the correlation observed in the first regime.

The net volume of the gulley and fan (volume accumulated in the fan minus volume eroded in the gulley) grew to 5.24e4 m 3 from January to June 2019, and then fell slightly to 4.0e4 m 3 by August, although errors are high and this may simply indicate that rockfalls had ceased. As the integration area was constant, the error in the volume was also constant at 1.20e4 m 3 , however the estimate of the uncertainty is only based on the standard deviation of four 'zero volume' gulley and fan dDEMs, so is itself likely quite uncertain. We measured the net volume of the gulley and fan in the hope that the erosion and deposition of older material cancels out leaving the net contribution of fresh volcanic material, however it seems likely that redeposited older material will not compact back down, and so there will be a net positive volume even before new volcanic material comes into consideration. The uncertain contribution of snow to these volumes also makes them suspect. However we include them to show that, in principle, the gulley and fan have a net positive volume that could account for a substantial part of the apparent decline in dome growth in June 2019 -i.e., rather than continuing to contribute to dome growth, it is at plausible that extrusion of further material triggered collapse and diverted material out of the crater instead during this period. geometry constrained by a conical-shaped crater and used to model the time evolution of the distance between the dome and crater edge (Fig. 7). TADR is a primordial parameter to monitor the evolution of effusive eruptions and yet is not always an easy parameter to determine with sufficient accuracy, precision, or temporal resolution. TADR estimation combined with prior knowledge of a volcano allows hazard assessments based on semi-quantitative event trees [START_REF] Newhall | Chapter 8 -Using Multiple Data Sets to Populate Probabilistic Volcanic Event Trees[END_REF][START_REF] Ogburn | The association of lava dome growth with major explosive activity (VEI ≥ 4): DomeHaz, a global dataset[END_REF]. Eruption parameters such as extrusion rate and magma composition integrated with

(1) a correct statistical analysis of global dome growth episodes database [START_REF] Newhall | Explosive activity associated with the growth of volcanic domes[END_REF][START_REF] Ogburn | DomeHaz: Dome-forming eruptions database[END_REF]; (2) the recent eruptive record and geological background of the volcano; (3) an adequate interpretation of instrumental and remote sensing monitoring and surface activity; and (4) the use of numerical simulation for modelling volcanic processes, allows for a trusted hazard assessment to be reached in complex risk scenarios situations and these were the parameters used by SERNAGEOMIN's frequent reports on the volcanic activity to forecast future possible scenarios (http://sitiohistorico.sernageomin.cl/volcan.php?iId=32). At Nevados de Chillán the complexity stems from the fact that several touristic facilities are exposed to volcanic risks, the high seasonal population is unfamiliar with volcanic hazards, few roads are available in case of evacuation and the eruption period is very long-lived (started in December 2015).

The rate of dome growth at Nevados de Chillán is one of the lowest ever reported. Table 2 shows the average growth rate at other well studied lava domes. By comparison, the maximum growth rate at Nevados de Chillán of ~1500 m 3 /day is about two orders of magnitude slower than most. Comparison with the global dome growth episodes database of [START_REF] Ogburn | DomeHaz: Dome-forming eruptions database[END_REF], reveals that only 12 out of 96 events (12,5%) with calculated time average extrusion rate have comparable values for TADR (<0.1 m 3 /s). Nevertheless, comparing with historical eruption at Nevados de Chillán, at the same subcomplex, the current eruptive cycle P a g e | 34 exhibits similar TADR with the last two historical magmatic eruption. The 2008 Sebastian lava flow, a viscous and dacitic lava flow erupted over an 8 month period with a mean TADR of 0.075 m 3 /s and no associated explosive activity [START_REF] Coppola | The 2008 "silent" eruption of Nevados de Chillán (Chile) detected from space: Effusive rates and trends from the MIROVA system[END_REF](Coppola et al., ) and the 1976(Coppola et al., -1986 Arrau eruption with a mean TADR of 0.005 m 3 /s [START_REF] Naranjo | Nuevos antecedentes sobre la evolución Cuaternaria del Complejo Volcánico Nevados de Chillán[END_REF]. As seen from the data in Table 2, there are no correlations between the growth rates and final dome volumes, so the lava dome growth rates are not informative of the final dome volumes. The time evolution of the dome growth rate observed at Nevados de Chillán is broadly similar to the one observed at other volcanoes. During the 2004 -2005 eruption of Mount St. Helens for instance the dome growth rate dropped from about 9 to 3 m 3 /s within the first three months [START_REF] Schilling | Use of digital aerophotogrammetry to determine rates of lava dome growth[END_REF]. Similar patterns where observed in the growth rates of the 2009 [START_REF] Diefenbach | Photogrammetric monitoring of lava dome growth during the 2009 eruption of Redoubt Volcano[END_REF][START_REF] Miller | Dome growth and destruction during the 1989-1990 eruption of redoubt volcano[END_REF][START_REF] Miller | Dome growth and destruction during the 1989-1990 eruption of redoubt volcano[END_REF][START_REF] Miller | Dome growth and destruction during the 1989-1990 eruption of redoubt volcano[END_REF]) lava domes at Redoubt volcanoes, 1991 lava dome at

Unzen [START_REF] Jaupart | Exogenous and endogenous growth of the Unzen lava dome examined by satellite infrared image analysis[END_REF][START_REF] Nakada | Endogenous growth of dacite dome at Unzen volcano (Japan), 1993-1994[END_REF] and 1980 lava dome at Mount St. Helens [START_REF] Fink | Intrusive and extrusive growth of the Mount St Helens lava dome[END_REF][START_REF] Swanson | Regularities in Growth of the Mount St. Helens Dacite Dome[END_REF].

The 2017-2019 eruptive episode at Nevados de Chillán is a striking example of effusive and explosive activities taking place simultaneous at a silicic volcano. The inception of dome formation and growth at the surface was accompanied by a sharp increase in explosive activity, averaging less than 5 to 10 eruptions per day prior to December 2017 to reach more than a hundred explosions in a day on 09 January 2018. Eruption of silicic volcanoes commonly show transitions from episodes of explosive activity to regimes of dome formation and vice versa (e.g., [START_REF] Fink | The Emplacement of Silicic Domes and Lava Flows[END_REF][START_REF] Nguyen | Explosive to effusive transition during the largest volcanic eruption of the 20th century (Novarupta 1912, Alaska)[END_REF]. The transition from one style of activity to the other is thought to occur due to changes in magma ascent rate (e.g., [START_REF] Dingwell | Volcanic Dilemma: Flow or Blow?[END_REF][START_REF] Martel | Timescales of bubble coalescence, outgassing, and foam collapse in decompressed rhyolitic melts[END_REF][START_REF] Woods | Transitions between explosive and effusive eruptions of silicic magmas[END_REF] or gas loss through permeable conduit walls (e.g., [START_REF] Eichelberger | Non-explosive silicic volcanism[END_REF][START_REF] Jaupart | Exogenous and endogenous growth of the Unzen lava dome examined by satellite infrared image analysis[END_REF]. Neither of these processes however easily explains how concurrent effusive and explosive activity can P a g e | 35

be taking place at a single vent. We envision that the conduit at Nevados de Chillán must be filled by a magma penetrated by an extensive fracture network allowing overpressure to be periodically released by ash venting or explosions accompanying magma ascent and dome growth at the surface (Fig. 9). Both processes are operating at different timescales but appear to be strongly positively correlated; periods of heightened explosion frequency corresponding to the periods of heightened dome growth rate and vice versa.

The strong first-order positive correlation between effusion rate and the frequency of explosions, long-period earthquakes and tremors (Fig. 5) suggests (1) that these metrics all tract the same physical process, presumably reflecting variations in overpressure of the magmatic system and (2) that one metric may be used to estimate another. As shown in this study and previous ones, measuring lava effusion rate in real time is not an easy endeavour for any volcano observatory. The strong correlation observed especially with the frequency of long-period events (R 2 =0.9) suggests that LP event frequency might be used as a simple way of estimating effusion rate in real time as the eruption continues. A similar observation was made after the 2018 Kīlauea eruption were a positive correlation (R 2 =0.6) between effusion rate and real-time seismic amplitude was observed [START_REF] Patrick | Cyclic lava effusion during the 2018 eruption of Kīlauea Volcano[END_REF]. In details however, the relationship between effusion rate and seismic activity at Nevados de Chillán is not a simple and continuous one and three regime can be distinguished showing strong increase in seismic events with increasing effusion rate at low effusion rates (≤ 100 m 3 /day), no noticeable changes in seismic activity for effusion rates between ~100 and ~900 m 3 /day, and again an increase in the number of seismic events with effusion rate at higher effusion rates (~900 to ~1700 m 3 /day) but with a different slope from the first regime (Fig. 5). Following the observation period reported in this manuscript the effusion rate increased again with a lava flow originating from the lava dome flowing down the northern flank. Then again, a strong positive correlation between effusion rate and the number of explosions, tremors and longperiod events was noticed but again with a different slope. The question we are left with is what physical process are these different regimes a reflection of? Regime 3, with effusion rates between ~900 to ~1700 m 3 /day, corresponds to the initiation of dome growth and includes the months from January to May 2018 (Fig. 4), regime 2 with effusion rates between ~100 to ~900 m 3 /day correspond to the period from June 2018 to April 2019 when effusion rate was decreasing very slowly (Fig. 4) with seismic activity remaining constant while regime 1 corresponds to the months of May 2019 to July 2019 when the effusion rate decreased faster and the number seismic events also dropped (Fig. 4). As the eruption continues and variations in effusion rates are recorded, we will be able to determine if these regimes are replicated or if a strong hysteresis emerges.

The disconnect in terms of ascent processes between the gas phase (sporadically escaping via a fracture network) and the magma might partly explain the extremely slow dome growth rate, much of the overpressure being accommodated by the gas escape leaving little driving force for the magma to ascend. Assuming a conduit radius of 5 meters, our dome growth rates can be translated to ascent rate at the top of the conduit (base of the dome) ranging from 20 meters per day in early 2018 to 9 meters per day in June 2019 (corresponding to decompression rates in the order of 10e-6 Mpa/s). Considering a smaller conduit would lead to higher ascent rate estimates (up to ~130 m/day for a 2 m radius conduit) which would remain much slower than most other well studied volcanic eruptions (see Fig. 9 in [START_REF] Shea | Bubble nucleation in magmas: A dominantly heterogeneous process?[END_REF], Fig. 5 in Cassidy et al., 2018and Fig. 13 in Moussallam et al., 2019 for compilations of ascent rates estimates).

Historical activity records at Nevados de Chillán include eruptive cycles with durations of years to decade [START_REF] Naranjo | Nuevos antecedentes sobre la evolución Cuaternaria del Complejo Volcánico Nevados de Chillán[END_REF]. This long-lived eruptive behaviour has been seen at P a g e | 37

other volcanoes with lava dome eruption such as Popocatépetl [START_REF] Mendoza-Rosas | Statistical analysis of the sustained lava dome emplacement and destruction processes at Popocatépetl volcano, Central México[END_REF] (also see global compilation by [START_REF] Wolpert | The longevity of lava dome eruptions[END_REF]. In addition, Holocene activity at Volcán Viejo (Las Termas Subcomplex) includes Vulcanian and sub Plinian episodes with 4 fall pumice and at least 3 pyroclastic density currents deposits recognized [START_REF] Dixon | The geology of Nevados de Chillán volcano, Chile[END_REF].

The wide range of past eruptive styles at Nevados de Chillan hence suggest that our two "worst-case" scenarios modelled here to map the maximum extent of PDC are realistic scenarios. The ongoing mostly effusive eruption could be but a stage in a bigger cycle, with much more explosive events remaining a future possibility. Monitoring the TADR and seismic activity at Nevados de Chillán in near-real time hence remain a priority to forecast and prepare for future activity. Further petrological studies would also help in this endeavour. Table 2: Growth rates of lava domes, data are from the compilation of [START_REF] Yokoyama | Growth rates of lava domes with respect to viscosity of magmas[END_REF] with 607 additional data from [START_REF] Miller | Dome growth and destruction during the 1989-1990 eruption of redoubt volcano[END_REF], ** [START_REF] Diefenbach | Photogrammetric monitoring of lava dome growth during the 2009 eruption of Redoubt Volcano[END_REF], and *** Schilling et al., 608 (2008) We have used aerial photogrammetry and Pleiades-1 tri-stereo imagery to produce ultra-highresolution DEMs tracking the growth of a lava dome at the Nevados de Chillán volcanic complex. We found that (i) The dome growth started as fairly stable over the first three months before decreasing over time (ii) The growth rate is orders of magnitude slower than most other lava domes (iii) Prediction of the time at which the lava dome should reach the crater rim proved accurate (iv) Modelling of dome-collapse-generated PDC showed no direct threat to populated areas but PDC generated by large explosions could travel towards populated area (v) The coexistence and simultaneity of effusive and explosive activity is uncommon and may contribute to explains the resulting exceptionally slow ascent rate (vi)

Positive correlation between dome growth rate, eruption frequency and seismic activity suggest that they all respond to the same driving force and that the frequency of LP event could be used as a proxy of effusion rate for operational monitoring.

Figure 1 :

 1 Figure 1: A: DEM of the Nevados de Chillán Volcanic Complex and surrounding valleys

  D o m e g r o w t h a t N e v a d o s d e C h i l l á n P a g e | 8

  D o m e g r o w t h a t N e v a d o s d e C h i l l á n P a g e | 10

  at given dates relative to our 21st Dec 2017 reference surface. For the gulley and fan, we took the 16th January 2019 DEM as our reference surface and differenced all successive DEMs with it. Histograms of the dDEMs over stable regions are shown in Fig S4 giving the Median Absolute Deviation (MAD) and standard deviation for each dDEM.

  D o m e g r o w t h a t N e v a d o s d e C h i l l á n P a g e | 11

  D o m e g r o w t h a t N e v a d o s d e C h i l l á n P a g e | 14

  D o m e g r o w t h a t N e v a d o s d e C h i l l á n P a g e | 15

Figure 2

 2 Figure 2 (panel 1 of 2): DEMs of the summit area after helicopter observation flights and

Figure 2 Figure 3 :

 23 Figure 2 cont (panel 2 of 2)

  2018 to January 2019, followed by another period of decreasing activity from January 2019 to June 2019 and increasing again in July to August 2019. The number of explosions at the crater per day (Fig. 4C) also follows the same temporal evolution, averaging 60 to 40 explosions per day in the first three months before stabilising around 35 explosions per day from April 2018 to April 2019, dropping below 20 eruptions per day in June 2019 before increasing sharply in July-August 2019. A direct comparison between our modelled mean monthly effusion rate and measured mean monthly numbers of seismic events is shown in P a g e | 21

D

  o m e g r o w t h a t N e v a d o s d e C h i l l á n P a g e | 22

Figure 4 :

 4 Figure 4: A. growth of dome and gulley-and-fan; red points are dome volumes, green points net volume of gulley and fan, blue tint is posterior

Figure 5 :

 5 Figure 5: Mean monthly modelled effusion rate at Nevados de Chillan for the period of

Figure 6 :

 6 Figure 6: A: Time evolution of the minimum horizontal distance between the foot of the lava

Figure 7 :Figure 8 :

 78 Figure 7: Time evolution of the minimum horizontal distance between the foot of the lava

D

  o m e g r o w t h a t N e v a d o s d e C h i l l á n P a g e | 36

DFigure 9 :

 9 Figure 9: Schematic depiction of the magmatic conduit at Nevados de Chillán, perforated by

Table 1 :

 1 Volume, Time Averaged Discharge Rates and Effusion Rates for the dome and gulley DEM timeseries. The uncertainties for a given value are given to the last two significant figures using concise notation (i.e. 5.04(13)e+04 = 5.04e4 ± 0.13e4). TADR values are given for the date at the end of the interval they are measured over. Effusion rates are instantaneous discharge rates predicted by the GP model for the specified date. The DEMs used as datums for estimating volumes of the dome and gulley are tinted red and green respectively.

	Date	Type	Volume m 3	TADR m 3 day -1	Effusion Rate	Gulley volume m 1	Gulley TADR
					m 3 day -1		m 3 day -1
	2017-12-21	Aerial	0.0(1.0)e+03	-	-	-	-
	2018-01-09	Aerial	2.49(17)e+04	1.31(10)e+03	1.73(11)e+03	-	-
	2018-01-12	Aerial	2.98(18)e+04	1.61(82)e+03	1.72(10)e+03	-	-
	2018-01-16	Aerial	3.57(19)e+04	1.48(64)e+03	1.694(88)e+03	-	-
	2018-01-23	Aerial	5.04(23)e+04	2.11(43)e+03	1.652(71)e+03	-	-
	2018-02-07	Aerial	6.86(27)e+04	1.21(24)e+03	1.552(44)e+03	-	-

  . 609

	610	Volcano Soufrière St. Vincent Colima Lamington Redoubt Redoubt* 1989-1990 Year 1979 1998 1951 1987 Redoubt** 2009 Soufrière Hills 1996 Turumai 1909 Bezymianny 1956 Shiveluch 1980 Pelée 1902 Mt. St. Helens 1980 Mt. St. Helens*** 2004-2005 Santiaguito 1922 Popocatépetl 1996 Unzen 1992 Usu 1944 Novarupta 1912 5. CONCLUSIONS P a g e | 40 Final lava dome volume (x10 7 m 3 ) Average growth rate (m 3 /day) 3.5 5.37E+05 0.04 3.80E+05 6 1.58E+06 6.8 2.00E+05 8.8 5.01E+05 7.2 7.94E+05 6.8 2.00E+05 1.5 3.80E+06 4.2 1.41E+06 1 1.86E+05 3.7 3.16E+05 1 8.32E+05 7.3 1.58E+05 20 1.58E+05 1.1 1.82E+05 5 5.01E+05 4.4 1.12E+05 0.5 2.51E+04

ACKNOWLEDGEMENTS

Pleiades images were acquired as part of the CNES' ISIS scientific support program (proposal ID: 2018-012-Sci, PI: YM), © CNES (2018-2019) and Airbus DS (2018-2019), all rights reserved, commercial use prohibited. This work was supported by public funds received in the framework of DINAMIS-GEOSUD, a project (ANR-10-EQPX-20) of the program "Investissements d'Avenir" managed by the French National Research Agency. We are very gratefully to two anonymous reviewers for their constructive comments and to Jean-Philippe Avouac for editorial handling.