David Castro-Cruz

Julie Régnier
email: julie.regnier@cerema.fr

Etienne Bertrand
email: etienne.bertrand@cerema.fr

Françoise Courboulex
email: courboulex@geoazur.unice.fr

A new parameter to empirically describe and predict the non-linear seismic response of sites derived from the analysis of Kik-Net database

Keywords: Non-linear site behavior, site effect, Kik-net, Ground motion prediction, Kumamoto Earthquake Highlights

We evaluate the non-linear site response at the stations of the Japanese Kik-Net network by computing the ratio between the Fourier spectra of the recordings at the surface and at the downhole station. When the amplitude of the input signal increases, we observe a shift of the resonance peaks towards lower frequencies, characteristic of non-linear soil behavior. We propose a new parameter called fsp to characterize this shift. We observe that fsp is a good proxy for the shear modulus reduction by comparing fsp and the modulus reduction curve obtain by both laboratory test measurements and numerical approaches. The fsp is used to correct the linear site response and surface ground motion in frequency domain in order to take into account the non-linear soil behavior. We show that the procedure we established reduces the error in predicting the main frequency peak by 50% compared to an elastic evaluation of the site response during the 2016 Kumamoto earthquake (Mw=7.1).

Introduction

It is widely recognized that soft soils can amplify dramatically the seismic motion compared to rock reference site because of seismic waves being trapped in sub-surface sedimentary layers [START_REF] Bard | The two-dimensional resonance of sediment-filled valleys[END_REF][START_REF] Cranswick | Site response across downtown Santa Cruz, California[END_REF][START_REF] Lermo | Site effect evaluation using spectral ratios with only one station[END_REF]. The seismic site response depends on the site geometric configuration, the nature of the soil and the incoming ground motion [START_REF] Kramer | Geotechnical earthquake engineering[END_REF]. For weak ground motion, the combination of a complex site configuration and various incident wave field induces a variability on the site response [START_REF] Thompson | Impediments to Predicting Site Response: Seismic Property Estimation and Modeling Simplifications[END_REF]. For strongest ground motions, the soil non-linear behavior can, in addition, have a strong influence on the site response [START_REF] Rubinstein | Nonlinear Site Response in Medium Magnitude Earthquakes near Parkfield, CaliforniaNonlinear Site Response in Medium Magnitude Earthquakes near Parkfield, California[END_REF][START_REF] Régnier | Assessing Nonlinear Behavior of Soils in Seismic Site Response: Statistical Analysis on KiK-net Strong-Motion Data[END_REF]. An accurate prediction of strong ground motion on sedimentary site will, therefore, require the consideration of non-linear soil behavior.

In laboratories, under cyclic and seismic loading, the soil non-linear behavior is illustrated in terms of hysteretic loops in the stress-strain plan [START_REF] Zeghal | Analysis of Site Liquefaction Using Earthquake Records[END_REF][START_REF] Assimaki | Quantifying Nonlinearity Susceptibility via Site-Response Modeling Uncertainty at Three Sites in the Los Angeles Basin[END_REF]. The effects of non-linearity in the soil mechanical behavior are generally a reduction of the shear modulus with increasing strain and an increase of the dissipation of energy within the material. Consequently, the hysteretic behavior of the soil can be recovered knowing the shear modulus reduction curve (G/Gmax) and the damping (ξ) curves .

For the transfer function, that quantify the site response, the effects of the non-linear soil behavior leads in a shift of the resonance frequencies of the site towards lower frequency, induced by a decrease of the apparent shear wave velocity, associated in most of the cases to an attenuation of the high frequency amplifications [START_REF] Régnier | Empirical Quantification of the Impact of Nonlinear Soil Behavior on Site Response[END_REF]. Nevertheless, an increase of the amplification at high frequency has been also sometimes reported, which was associated with a decrease of the pore water pressure when cyclic mobility occurs [START_REF] Bonilla | Hysteretic and Dilatant Behavior of Cohesionless Soils and Their Effects on Nonlinear Site Response: Field Data Observations and Modeling[END_REF] or soil hardening [START_REF] Pavlenko | Possible Mechanisms for Generation of Anomalously High PGA During the 2011 Tohoku Earthquake[END_REF]. Many studies show that these differences are larger when the intensity of the incident motion increases (e.g. [START_REF] Aguirre | Nonlinearity, liquefaction, and velocity variation of soft soil layers in Port Island, Kobe, during the Hyogo-ken Nanbu earthquake[END_REF][START_REF] Iai | Response of a Dense Sand Deposit during 1993 Kushiro-Oki Earthquake[END_REF][START_REF] Wen | Non-linear soil response in ground motions[END_REF]).

The non-linear soil parameters can be characterized by various methods from laboratory testing on soil samples to lower cost method based on correlations between non-linear parameters and other soil parameters (e.g. [START_REF] Ishibashi | Unified dynamic shear moduli and damping ratios of sand and clay[END_REF][START_REF] Vucetic Mladen | Cyclic Threshold Shear Strains in Soils[END_REF]). In-situ measurements can also be used in addition to laboratory tests to infer these non-linear parameters. Some authors use accelerometric data to determine the shear modulus reduction and damping curves. These parameters are often estimated by finding the best numerical model that reproduces the accelerometric time histories. It has been applied mainly in Taiwan and in Japan (e.g. [START_REF] Zeghal | Lotung Downhole Array. II: Evaluation of Soil Nonlinear Properties[END_REF][START_REF] Glaser | System identification estimation of soil properties at the Lotung site -ScienceDirect[END_REF][START_REF] Pavlenko | Estimation of Nonlinear Time-dependent Soil Behavior in Strong Ground Motion Based on Vertical Array Data[END_REF]). Other authors use interferometry between recordings at different depth to derive the instantaneous wave propagation velocity in the media depending on the input motion intensity [START_REF] Nakata | Estimating near-surface shear wave velocities in Japan by applying seismic interferometry to KiK-net data[END_REF][START_REF] Bonilla | Prediction of non-linear site response using downhole array data and numerical modeling: The Belleplaine (Guadeloupe) case study[END_REF].

Another way to study soil non-linear behavior is to compare site-response curves computed from weak motion and strong motion. In Noguchi and Sasatani [START_REF] Noguchi | Quantification of degree of nonlinear site response[END_REF] the ratio between the Fourier spectrum of records at the surface and at downhole (BSR) are computed. They compared BSR from weak and strong ground motion by the summation of the differences at each frequency between each BSR, parameter that was called DNL (Degree of Non-linearity). They showed a link between the nonlinear effects (defined as DNL), and an intensity parameter of the shaking (PGA). In Lussou et al. [START_REF] Lussou | Quantification of soil non-linearity based on simulation[END_REF] and Régnier et al. [START_REF] Régnier | Assessing Nonlinear Behavior of Soils in Seismic Site Response: Statistical Analysis on KiK-net Strong-Motion Data[END_REF], several parameters are proposed to quantify those differences, such as an estimation of the modification of the BSR curves and a measurement of the frequency shift. In Field et al. [START_REF] Field | Nonlinear ground-motion amplification by sediments during the 1994 Northridge earthquake[END_REF] ratios between spectral ratios were computed using synthetic rock reference seismograms. Similarly, in Régnier et al. [START_REF] Régnier | Empirical Quantification of the Impact of Nonlinear Soil Behavior on Site Response[END_REF] the ratio of the non-linear to linear spectral ratios are analyzed, but using vertical arrays. In all cases, the authors found that the non-linear to linear site response discrepancy is related to the frequency range and increase with increasing PGA of the incident ground motion.

Following a similar approach, we compare weak motion and strong motion site responses in order to propose a methodology for correcting the linear site response in order to reproduce the effects of non-linear soil behavior based on the analysis of recordings from vertical arrays of the Japanese Kiban Kyoshin network (KiK-Net). After a description of our approach, an application will be detailed on the recordings of the 2016 Kumamoto Japanese earthquake. Subsequently, for all of the selected seismograms a signal processing based on Boore [START_REF] Boore | On Pads and Filters: Processing Strong-Motion Data[END_REF] recommendations was applied. It consists in: (1) removing the mean, (2) removing the data before the first zero-crossing point, (3) Applying an Hanning's window to improve the trimming process of the signal [START_REF] Oppenheim | Discrete-time signal processing[END_REF], and finally (4) Applying two times a high pass filter of 2nd order with fc=0.1 Hz. This processing allows to remove the low-frequency noise and to ensure compatibility between acceleration, velocity and displacement time series and Fourier and response spectra. In particular, effects of spurious low frequency can be removed from the data through this procedure.

The empirical site response is assessed on the horizontal components by computing the Borehole Spectral Ratios (BSR), defined as the ratio between the geometric average of the Fourier amplitude spectra of the horizontal components at the surface and at downhole (Equation ( 1)).

( ) = 2 + 2 ℎ 2 + ℎ 2 (1)
Where EW and NS are the amplitude of the Fourier transform of the accelerograms for East-West, and Northsouth horizontal components respectively. In this study, we applied a Konno-Ohmachi smoothing [START_REF] Konno | Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor[END_REF] to each spectrum, improving the comparison between spectra and avoiding zero values at the denominator of (Equation ( 1)). BSR shows, in the frequency domain, how the seismic signal varies between the downhole sensor and the one at the surface. BSR curves indicate the site resonance frequencies f (n) if the downhole sensor is located at the main sediment to substratum interface. Otherwise, the curve includes pseudoresonances that are induced by the destructive interferences of the down-going waves with the incident ones at the downhole sensor [START_REF] Field | A comparison and test of various site-response estimation techniques, including three that are not reference-site dependent[END_REF]. However, the comparison of weak and strong site response remains relevant to characterize non-linear site behavior.

Method

Our method defines a parameter called fsp that characterize the non-linear site behavior. To illustrate our methodology, we detail the whole procedure on one site of the KiK-net database (IBRH11) that has recorded both weak and strong ground motions. Then, we compare the parameter fsp to laboratory data (shear modulus reduction curves) obtained from cyclic tri-axial test at one Kik-net site (KSRH10), that was characterized during a previous PRENOLIN project [START_REF] Régnier | PRENOLIN: International Benchmark on 1D Nonlinear Site-Response Analysis-Validation Phase Exercise[END_REF]. Finally, the parameter is tested analytically by using non-linear site response calculations on simple site configuration (1-D, non-linear, mono-layer cases).

fsp: a new parameter to characterize non-linear site behavior

In order to explain our methodology, we detail the whole method on one site of the KiK-net database that has recorded a large panel of earthquakes from weak to strong motions. Then, we will generalize the procedure to all of the stations. Station IBRH11 has recorded 515 earthquakes and it is located in the prefecture of Ibaraki Figure 3 illustrates four BSR computed at IBRH11 from ground motions with different levels (PGAdownhole). We observe that weak ground motions share a similar BSR. This is because for weak ground motion the site response is linear [START_REF] Aguirre | Nonlinearity, liquefaction, and velocity variation of soft soil layers in Port Island, Kobe, during the Hyogo-ken Nanbu earthquake[END_REF]. The variability that can still be observed in the site response for these small seismic solicitations is mainly due to complex site geometry associated to the distribution of earthquake sources around the site [START_REF] Thompson | Impediments to Predicting Site Response: Seismic Property Estimation and Modeling Simplifications[END_REF]. For stronger ground motions, the peak frequencies occur at a lower frequency bandwidth. We interpret this shift as a direct effect of the shear modulus reduction of the soil when the seismic solicitation increases. Furthermore, the decrease in the amplitude of BSR is another effect of non-linear soil behavior linked to the increase of damping with shear strain, but in the following we will only discuss about the frequency shift impact. The weak ground motions for which the site behaves linearly are selected based on their maximal peak accelerations. Indeed, the ground motions with a PGAdownhole from 10 -4 m/s 2 to 6. 10 -3 m/s 2 are considered. Figure 4 presents the BSRlinear with a dotted line. It is the arithmetic average of all of the weak ground motions BSR (grey curves). The black line represents the BSR computed from a strong ground motion with PGAdownhole of 2.6 m/s 2 . It is easy to appreciate that the resonance peaks of the BSR from this strong ground motion are shifted to lower frequencies and that their amplitude are the smallest. BSRlinear is compared to each BSR from all of the recorded ground motions, allowing us to define a parameter that characterizes the observed frequency shift. The logarithmic frequency shift is the gap in logarithmic scale between BSRlinear and BSR. In linear scale, it is a coefficient that changes the frequency scale. The algorithm to find this logarithmic shift minimizes the misfit between BSRlinear and BSR as defined in the Equation (2). Note that the misfit is weighted by the logarithmic sampling.

= ! " # $%&' ( ̅ * , - . ̅ /" ∆1 ∆1 = 234 5 ( 65 / ) 78 ̅ = 0.5 • ( 65 + ) (2)
In the Equation ( 2) Ls is the logarithmic shift applied to BSRlinear. The misfit is defined as a discrete approximation of the area between the shifted BSRlinear and BSR, considering a logarithmic scale as the length of the base (∆x). The computation is done over a frequency window going from 0.3 Hz to 30 Hz.

Finally, we define a frequency shift parameter, so called fsp, as the square of the Ls, which produces the minimum value of misfit (fsp=Ls 2 when misfit is minimized). The fsp is a coefficient that is applied to the linear resonance frequency to obtain the non-linear ones for a specific ground motion. It means that if no shift is needed to fit both curves, both Ls and fsp are equal to one. If BSRlinear needs to be logarithmically shifted to higher frequencies to fit BSR, fsp will be higher than one. In the opposite case, fsp will be lower than one. Nonlinear soil behavior is expected to shift the BSR to lower frequency range and therefore linked to an fsp below one.

fsp is computed for every recording collected at station IBRH11. The results (Figure 5) shows the evolution of fsp in function of the ground motion intensity, expressed in terms of the PGA recorded at the downhole station. For small PGA, fsp is close to 1, but as the solicitation level increases, the value starts to decrease. The decrease is gentle for PGA smaller than 0.1 m/s 2 but is growing rapidly for larger seismic solicitations. A hyperbolic curve is used to fit the fsp values (Equation ( 3)), with a formulation equivalent to the one used to describe the modulus reduction curves [START_REF] Duncan | Nonlinear analysis of stress and strain in soils[END_REF][START_REF] Ishihara | Soil Behaviour in Earthquake Geotechnics[END_REF].

= 5 56 <=> ?@ABC@DE <=> FEG ?@ABC@DE (3)
Where HIJ K%L MNO$PN#% is a parameter that describes the hyperbolic function for each site. It is equal to the PGAdownhole corresponding to fsp=0.5, this parameter being the counterpart of γref in the hyperbolic models and the only parameter that is needed to fully describe the hyperbolic function. In order to quantify the variability of fsp, a single standard deviation is computed for each station considering the whole number of data. The dashed lines on Figure 5 represents the mean plus or minus one standard deviation.

Here, PGAdownhole can be considered as an estimator of the strain level that the site can experience during each earthquake. Others proxy of strains have already been proposed in the literature. For example, Chandra [START_REF] Chandra | PGA-PGV/Vs considered as a stress-strain proxy for predicting nonlinear soil response[END_REF] proposed PGV/VS30 as a strain proxy. We decided to use first the PGA recorded at the downhole station instead because the integration from acceleration to velocity may add some uncertainties in the PGV calculation. Also PGA is a relevant parameter for describing the amount of non-linearity a soil may produce [START_REF] Régnier | Assessing Nonlinear Behavior of Soils in Seismic Site Response: Statistical Analysis on KiK-net Strong-Motion Data[END_REF].

fsp parameter compared to the shear modulus reduction curves

We compute the fsp parameter for the earthquakes recorded at station KSRH10. This site is situated in the Hokkaido prefecture, has recorded 283 events, and is located on a D-site according to EC8 (Vs30 of 212 m/s). This station is particularly interesting because it has been well characterized by in-situ and laboratory data measurements performed during the PRENOLIN project [START_REF] Régnier | PRENOLIN: International Benchmark on 1D Nonlinear Site-Response Analysis-Validation Phase Exercise[END_REF]. We plot the fsp values computed for this site against the PGV at the surface divided by Vs30, as a proxy for strain, to compare with the shear modulus reduction curves obtain from cyclic tri-axial tests performed at this station (Figure 5). We observe that the fsp values are in the range of the shear modulus reduction curves. From this comparison we are convinced that fsp is a good proxy for the shear modulus reduction with shear strain.

Comparison with analytical and numerical modeling

We also test fsp analytically and with numerical modeling. First, we analyze the fsp analytically in a very simple case and then we apply our calculations to a more complex wave propagation in truly non-linear soil model.

In linear domain with 1-D configuration of a mono-layer soil and considering vertical incidence of the shear waves, the resonance frequencies (f {n} ) depend on the shear velocity (Vs) and the layer thickness (H) through f {n} =Vs(1+2n)/(4H). The shear velocity (Vs) can be expressed as a function of the shear modulus as Q R = SI T ⁄ where G and ρ are the shear modulus and the density respectively. The shear modulus of the layer can thus be derived knowing the associated resonance frequencies using the Equation (4).

I = ( 2 • V • {$} (0.5 + 8) Y , Z T (4)
In the equivalent linear approximation (EQL), the nonlinearity of the soil behavior is taken into account in an iterative process that adjusts the elastic properties to the level of strain induced in the layer, knowing both the modulus reduction and the attenuation curves. In this frame, the shear modulus at a linear strain rate (Gmax) can be compared to the shear modulus at larger strain (G), by the ratio between both. Using the Equation (4) for both modulus, the ratio will be expressed as followed:

I I [&\ = ] {$} # $%&' {$} Y ^Z = (5)
The Equation [START_REF] Thompson | Impediments to Predicting Site Response: Seismic Property Estimation and Modeling Simplifications[END_REF] shows that the ratio of the shear modulus is proportional to the square root of the ratio between the linear resonance frequency flinear and the non-linear ones (whatever the order of the harmonic).

In logarithmic scale, this ratio (f/flinear) represents a logarithmic shift, and this shift can be found minimizing the Equation [START_REF] Cranswick | Site response across downtown Santa Cruz, California[END_REF]. Computing the square of the logarithmic shift that minimizes the Equation ( 2) is, once again, the definition of the fsp parameter. Therefore, in this very simple case, the shear modulus reduction is equal to the previously defined fsp.

In a fully non-linear model, it is not possible to derive an analytical formulation as for the equivalent linear approach. Therefore, the relationship between the shear module reduction with strain increase and the shift of the resonance frequencies towards lower frequencies must be analyzed using numerical simulation and numerical BSR calculation. To simulate the seismic response, we are using a 1-D fully non-linear approach implemented in CyberQuake software proposed by Modaressi et al. [START_REF] Modaressi | User's Man BRGM Fr[END_REF] and based on a plastic constitutive model with hardening based elastoplastic theory [START_REF] Hujeux | Une loi de comportement pour le chargement cyclique des sols[END_REF]. We compare the results to the ones computed with an equivalent linear model also implemented in CyberQuake. We are considering a model composed of a single layer of 40 m of thickness with a shear wave velocity of 200 m/s over an elastic bedrock. In the linear domain, the fundamental resonance frequency of this soil column is equal to 1.25 Hz. For the EQL, the non-linear properties are homogeneous along the depth and illustrated in Figure 7. For the fully non-linear approach, the soil mechanical parameters are also shown in Figure 7. Those parameters were chosen as they characterize a typical non-linear soil layer [START_REF] Lopez-Caballero | Nonlinear numerical method for earthquake site response analysis I -elastoplastic cyclic model and parameter identification strategy[END_REF]. The non-linear constitutive model implemented defines non-linear properties that depends on the confining pressure and consequently on the considered depth, while Vs remains the same in the whole layer. Figure 7 shows in green the modulus reduction curves and damping curves at different depths.

The bedrock is modeled as an infinite monotonic and elastic material with shear velocity of 1500 m/s, a density of 2000 kg/m 3 , and a Poisson coefficient of 0.3. A Gabor wavelet is used as input signal. This simple wavelet makes easier the evaluation of non-linear effects in the BSR. The Gabor's wavelet is defined with the Equation (6).

_( ) = J

`abcd e (fgfh)

i j b k3 (2l [ ( -3)) (6)
In the simulations, γ is taken as 3, to as 2.5s. The central frequency (fm=1.25 Hz) is chosen in agreement with the fundamental resonance frequency of the site and only frequencies lower than 2 Hz are analyzed, as the energy of the input signal at higher frequencies is too small to make the ratio computation relevant. A corresponds to the maximal amplitude of the signal. For the fully non-linear approach, the soil mechanical parameters are given in Table 1. They were chosen as they characterize a typical non-linear soil layer [START_REF] Lopez-Caballero | Nonlinear numerical method for earthquake site response analysis I -elastoplastic cyclic model and parameter identification strategy[END_REF]. In the linear domain, the fundamental resonance frequency of this soil column is equal to 1.25 Hz.

Table 1. Parameters of the soil layer for the 1-D numerical simulation.

Vs 1 (m/s) 200 Vp 2 (m/s) 374 ρ 3 (Kg/m 3 ) 1750 Φ 4 (°) 25 β 5 10 σci/σv 6 2 Ep 7 20 C 8 (kPa) 0 γelastic 9 1x10 -9 b 10 0.9 nr 11 0.4 ψ 12 (°) 20 1 Shear wave velocity (m/s) 2 Primary wave velocity (m/s) 3 Density (Kg/m 3 ) 4 Friction angle at totally mobilized plasticity (°) 5 Plastic modulus 6 Compaction ratio 7 Plastic stiffness coefficient 8 Cohesion (kPa) 9 Extent of the truly elastic domain 10 Shape parameter of the yield Surface 11 Numerical parameter of isotropic hardening 12 Slope of the characteristic line (°) αψ 13 0.8

We test different input signal amplitudes and calculate the BSR in each case. The maximum amplitudes of the signals are: 0.001, 0.01, 0.05, 0.1, 0.3, 0.7 and 1 m/s 2 , results are shown in Figure 8. For the smallest PGA, the computed BSR is equivalent to the analytical elastic case (black dash line). As expected for the larger input motions, the amplification peak is smaller, and it is shifted towards lower frequencies. The equivalent linear computation (Figure 8A) leads to resonance frequency shift larger than the fully non-linear approach (Figure 8B). Indeed, besides the limitations of EQL approaches at large strain [START_REF] Kaklamanos | Critical Parameters Affecting Bias and Variability in Site-Response Analyses Using KiK-net Downhole Array DataCritical Parameters Affecting Bias and Variability in Site-Response Analyses[END_REF], the non-linear soil properties differ in the two numerical approaches as mentioned above. Indeed, at the soil to bedrock interface, the non-linear soil properties are much more linear in the fully non-linear computation. In the linear equivalent approach, the effective strain (γeff) is used to calculate the equivalent shear modulus and damping at each iteration. It is defined from the maximal strain as γeff= 0.65γmax. In Figure 9a, the values of the G/Gmax ratio computed using γeff and γmax are compared to the appropriate fsp values. Taking into account γmax leads to a misfit between fsp and the ratio at high deformation level. Those differences are related to the EQL methodology who computes the non-linear parameters with γeff instead γmax. Considering G/Gmax at yeff, the frequency shift equals almost the shear modulus reduction. The small error can be attributed to the iteration scheme implemented in the linear equivalent method. The good fitting between G/Gmax and fsp curves shows that fsp parameter is able to quantify the amount of shear modulus reduction as equivalent linear method does in this very simple mono-layer case.

Figure 9b presents the G/Gmax ratio at different depths in function of the fsp parameter calculated from the BSR curves for a full non-linear approach. The G/Gmax ratio used in these graphs is computed considering the effective strain (γeff) at the middle of the layer for various depths in the non-linear computations (Figure 7).

From this figure, we observe that G/Gmax and fsp are not equivalent at all depths in the fully non-linear approach. Indeed, the curves showed on the graph depend on the depth considered for the computation of the shear modulus, as both the non-linear soil properties and the shear strain varies with depth. However, we notice that Equation ( 5) still can be applied to the curve computed at 25.5 m at depth. At this specific depth, the reason why fsp is not perfectly equal to the G/Gmax ratio on the whole fsp range could be related to the definition of γeff using an approximation of the shear strain during the shaking. At shallower depths, fsp overestimates the shear modulus decrease. On the contrary, at larger depths fsp underestimates slightly the G/Gmax ratio.

Prediction of non-linear surface ground motion using fsp curves

In order to go further, we aim to develop a methodology to predict the non-linear surface ground motion using the fsp curves previously computed. For a given site, knowing the PGA value at the downhole sensor, we propose to estimate fsp from the non-linear regression defined in the Equation ( 3) previously computed at the station. Combining fsp and BSRlinear, it is possible to estimate the BSR for the predicted ground motion. This estimation consists in applying the predicted frequency shift to BSRlinear as shown in Equation [START_REF] Régnier | Assessing Nonlinear Behavior of Soils in Seismic Site Response: Statistical Analysis on KiK-net Strong-Motion Data[END_REF].

m ( ) = 2 8 7 ( • S ) (7)
where m is the estimated borehole spectral ratio for the strong motion, BSRlinear the linear borehole spectral ratio, and fsp, the estimated frequency shift from the hyperbolic curve defined for the site (Equation ( 3)). PGArefdownhole is the parameter that quantifies the susceptibility of a site to show a non-linear behavior. The higher it is, the less non-linear the site is. We can observe that the variation between stations is significant. The PGArefdownhole go from 3.1x10 -1 m/s 2 for the station HRSH13 to infinite for AICH04. The hyperbolic trend previously describe is well suited for the fsp distribution in function of PGAdownhole at almost all of the stations. However, we can note that at station IWTH21, the hyperbolic model (blue line) is not necessarily the function that would produce the best fit.

In Table 2, we conclude that the higher Vs30 is, the smaller PGArefdownhole is, meaning that larger non-linear behavior is observed for site characterized by high Vs30. This has already been reported in previous papers (e.g. Régnier et al., 2014 [START_REF] Régnier | Assessing Nonlinear Behavior of Soils in Seismic Site Response: Statistical Analysis on KiK-net Strong-Motion Data[END_REF]) with the conclusion that high Vs30 sites, composed of a thin soft soil over a stiff shallow bedrock, create a very large impedance contrast that can induce strong amplification and large strain in the superficial layer during earthquakes.

4.2

Application to 2016 Kumamoto earthquake recordings, Mw 7.1

We apply our methodology to the Kumamoto Earthquake recordings. The earthquake mainshock occurred April 15th, 2016 in the south of Japan (island of Kyushu), with a magnitude Mw of 7.0. This earthquake initially started on a deep portion of the northern part of the Hinagu fault and then finished in the Futagawa fault [START_REF] Asano | Source rupture processes of the foreshock and mainshock in the 2016 Kumamoto earthquake sequence estimated from the kinematic waveform inversion of strong motion data[END_REF], it was the largest earthquake in Japan in 2016. In this section, the recorded ground motions at KIK-net sites (FKOH03, OITH11, and KMMH03) are compared to predicted ones using the methodology described previously. Figure 12 gives the location of the three sites in the Island of Kyushu, the location of the fault and the epicenter according to [START_REF] Yagi | Rupture process of the 2016 Kumamoto earthquake in relation to the thermal structure around Aso volcano[END_REF]. Station FKOH03 is located at 91 km of the fault located in Umi Fukuoka prefecture, with a Vs30 of 497 m/s (site EC8 type C). The bedrock is reached at 25m depth with a Vs of 2030 m/s and the downhole station is located at 100m depth. Station OITH11 is located at 72 km from the epicenter in Kokonoe Oita prefecture. It has been characterized with a Vs30 of 458 m/s (site EC8 type C) and the S wave velocity reaches 800 m/s starting from 65m depth, the downhole station is located at 160 m depth with a Vs of 1000m/s. Finally, KMMH03 is located in the prefecture of Kumamoto in Kikuchi at 27.7 km from the epicenter. Like the two other sites, it is a EC8 type C class with a Vs30 of 421 m/s. The bedrock was found at 80m with a Vs of 1200 m/s and the downhole station is located at 200m depth where Vs of reaches 2000m/s. We compute the fsp curves for the three sites (Figure 14) using 140, 57, and 111 earthquakes for station FKOH03, OITH11, and KMMH03 respectively leaving on the side the recordings of the main Kumamoto event.

On figure 14, we report together with the fsp curves, the PGA values at downhole recorded during the Kumamoto earthquake (orange vertical line). We note that the ground motions recorded during this event were larger than all the other events of the database for stations OITH11 and KMMH03. Therefore, for these stations, fsp is estimated by extrapolation of the fsp curves calculated on weaker motions. The performance of the extrapolation depends on the site and the available data. If no large or medium earthquakes have been recorded at a site, the non-linear regression is not well constrained at larger PGA and the estimated shift will be assorted with important uncertainties. On the contrary, at FKOH03, the fsp is assessed by interpolation leading to a better confidence. However, even here, the dispersion between the curve and the real data remains. The p derived from the BSRlinear with the predicted frequency shift are illustrated in Figure 15 in black thick lines along with the observed BSR in grey lines and the linear BSR in dotted line. For the three cases, the estimated p is closer to the observed BSR than to the linear BSR. This result shows that the shift correction we propose is relevant for large ground motions, when the nonlinear soil behavior is running. Furthermore, we observe that the shift in the BSR is well predicted even with an extrapolated fsp.

At station FKOH03 that experienced the weaker motion during this earthquake (PGAdownhole ~ 0.2 m/s 2 ), the correction on BSRlinear is very small and do not particularly change the frequency content of the motion. But at station OITH11, that recorded a larger PGAdownhole (about 2m/s 2 ), the correction is important leading in a prediction fiting very well the real surface data. On this station, the BSR has one main peak and a simple shape. Station KMMH03 also recorded a strong PGAdownhole but its BSR is more complex, with various peak frequencies (mainly two around 1Hz and 10Hz) over a large frequency band. This station is located within few kilometers of the active fault. Consequently, it is subject to near-source high frequency effects. The correction applied is very efficient on the lower frequency of the BSR but tends to shift a little too much the higher frequency part (> 10z). For this station, we observe that the shift shouldn't be the same for all frequencies. This could be due to others causes than the non-linear soil behavior alone. Table 3 illustrates the improvement offered by the proposed methodology in the predicting of strong ground motion. We compute the difference in frequency of the BSR peaks between the observation and the related prediction. From this table, it is clear that the differences are lower when the BSR is corrected by the fsp shift. Indeed, at OITH11 and KMMH03, the methodology improves the prediction of the frequency peak by around 50%. At FKOH03, the ground motion was not strong enough to show some non-linear soil behavior during the main event. However, even in this case, we can see that a little improvement is provided by our approach. The Fourier spectrum for the geometric average of the horizontal components is computed according to the Equation ( 1). The results are presented in Figure 16 for the three sites on which we are focusing here. The predictions of the surface ground motion are improved when the site response is corrected by the adequate frequency shift. Comparing the main peaks frequency position of the Fourier spectra at the surface, Table 4 shows the relative difference between observations and predictions. In the case of the site FKOH03, the difference is small since the target ground motion isn't strong enough to exhibit a non-linear behavior of the soil. But even in this case, the spectra obtained with the fsp correction is closer to the observed spectra.

For the two other sites, the prediction of the main peak frequency position is improved by more than 40% in both cases in comparison with the direct use of BSRlinear.

Table 4. Comparison of the main peaks in the spectra at the surface between the observation and the computed spectra using qrs p and BSRlinear.

FKOH03 OITH11 KMMH03

Asurface to p •Adownhole 4% 8% 0% Asurface to BSRlinear• Adownhole 10% 60% 44%

Discussion

To estimate precisely the ground motion at the surface in the time domain, an inverse Fourier transform should be applied to the estimation of the Fourier spectrum. Additionally, the same previously described process should be applied independently for each horizontal components of then ground motion. However, the phase modification due to the site effects should also be taken into account. This effect has been studied before, and usually, it is assumed that 1D site effect does not modify consequently the phase [START_REF] Fleur | Site Effects in Port-au-Prince (Haiti) from the Analysis of Spectral Ratio and Numerical Simulations[END_REF]. Since we are using borehole arrays configuration, the applicability of the same assumptions must be evaluated, and it will be part of future studies.

The application of the parameter fsp is very promising for the prediction of the non-linear ground motion. In a simple case, we demonstrate that the parameter fsp provides a good evaluation of the shear modulus reduction. In a full non-linear approach, the parameter provides an average quantification of the decrease of the shear modulus. For both cases, the parameter gives a good quantification of the impact of the non-linear soil behavior in the site response.

For now, the definition of the fsp parameter requires a vertical array of accelerometers and weak to moderate earthquake recordings to define the fsp curves. To overpass this limitation, we are currently working on the prediction of the fsp curves at non-instrumented sites but characterized by specific geotechnical parameters as Vs30 and f0 that has been demonstrated to be a good proxy parameter for describing the potential nonlinear site behavior of a site (e.g. [START_REF] Régnier | Empirical Quantification of the Impact of Nonlinear Soil Behavior on Site Response[END_REF] Régnier et al 2016). To predict the rock motion at the downhole station (with the effects of the down-going waves) several methods could be applied such as the Empirical Green's Function directly at downhole or with a correction method to transfer from surface to borehole as proposed in Cadet et al. [START_REF] Cadet | Site effect assessment using KiK-net data: Part 1. A simple correction procedure for surface/downhole spectral ratios[END_REF] or the use of numerical models.

Conclusions

We propose in this paper to quantify the effects of non-linear site behavior by measuring the frequency shift that it produces in the Borehole Spectral Ratio (BSR) from weak motion to stronger one. The logarithmic shift of BSR, that we called frequency shift parameter (fsp), is related to the intensity of the ground motion. For each site, we, therefore, defined fsp curves that relate the fsp parameter to the PGA at the downhole station following the classic hyperbolic model. We thus are able to estimate the non-linear effects of strong ground motions at sites using enough recorded data.

We showed that the frequency shift parameter (fsp) is a good proxy for non-linear site behavior, as the fsp curves could mimic the shear modulus reduction curves. This fitting is better accomplished when an appropriate parameter to approximate the shear strain in the soil column during the shaking is used.

For the earthquake of Kumamoto 2016, we estimated the non-linear site response and the Fourier spectrum at the surface, including the non-linear effects based on the linear site response and the corresponding rock downhole recording of the earthquake. The prediction provided very close results to the observations and, for sites with high nonlinearity, it improves in our case by more than 40% the frequency position of the main peak compared to a linear evaluation. The fsp curves provide, as shear modulus reduction curves do, an evaluation of the propensity of a site to develop non-linearity as well a possible way to quantify the level of non-linearity of the site under strong earthquakes.

Figure 1 .

 1 Figure 1. Magnitude/distance distribution of earthquakes recordings used in this study.

(Figure 2 )

 2 on a D-site according to the Eurocode 8 (EC8) with a Vs30 of 242 m/s ( = 30 ∑ ℎ / ) and a depth of 103 meters. The bedrock is reached at 30m depth and is characterized with a Vs of 2100 m/s. The downhole station is situated at 100m depth.

Figure 2 .

 2 Figure 2. Location and Vs profiles of stations IBRH11 and KSRH10 and the other stations (Japanese Kik-net network).

Figure 3 .

 3 Figure 3. BSR at IBRH11 station (KiK-net) for four earthquakes with different PGA at the downhole station.

Figure 4 .

 4 Figure 4. Definition of BSRlinear as the average of all of the weak ground motions (142 events) and comparison between BSRlinear with the BSR computed from one strong ground motion.

Figure 5 .

 5 Figure 5. fsp value against PGA at downhole (fsp curves). Black line best hyperbolic function that fits the data. Station IBRH11 -VS30: 242 m/s.

Figure 6 :

 6 Figure 6: Shear modulus reduction curves defined in the PRENOLIN project at station KSRH10 (green lines) at different depths from cyclic triaxial tests and fsp against PGV/ Vs30 calculated at the site. The color version of this figure is available only in the electronic edition.

Figure 7 .

 7 Figure 7. Modulus reduction curve and damping curve that are used in: Equivalent Linear approx. (black line-squares), and at different depths for the full non-linear model (green lines). The color version of this figure is available only in the electronic edition.

Figure 8 .

 8 Figure 8. BSR for the equivalent linear approach (a) and the fully non-linear simulation (b) for a different level of input solicitation. Black dash is the elastic response curve.

Figure 9 .

 9 Figure 9. Evolution of the shear decay reduction in function of fsp. Left: equivalent linear approach. Right: full non-linear simulation. Red dotted curve: G/Gmax = fsp. The color version of this figure is available only in the electronic edition.

Figure 11 .

 11 Figure 11. fsp curve for several stations that represents the general behavior of the data set. The blue line is the hyperbolic function that is fitted in each case. The color version of this figure is available only in the electronic edition.

Figure 12 .

 12 Figure 12. Island of Kyushu Japan with the location of the studied stations. The line and the star represent the location of the active fault and the epicenter for Kumamoto earthquake 15 th April 2016 [39].

Figure 13 .

 13 Figure 13. Vs profiles of evaluated sites for the earthquake of Kumamoto 2016. (a) site FKOH03. (b) site OITH11. (c) site KMMH03.

Figure 14 .

 14 Figure 14. fsp parameter in function of the PGAdownhole. (a) site FKOH03. (b) site OITH11. (c) site KMMH03. The vertical line represents the PGAdownhole generated by the Kumamoto 2016 earthquake at each site.

Figure 15 .

 15 Figure 15. Comparison between the BSRlinear for each site (black dashed line) computed on weak motion, the BSR computed on the Kumamoto earthquake recordings (orange line), and the correction of the borehole transfer function (qrs m ) that we propose (dark blue line). (a) site FKOH03. (b) site OITH11. (c) site KMMH03. The color version of this figure is available only in the electronic edition.

Figure 16 .

 16 Figure 16. Fourier spectrum at the surface of the horizontal components. Earthquake of Kumamoto 15 th April. (a) site FKOH03. (b) site OITH11. (c) site KMMH03. The color version of this figure is available only in the electronic edition.

Table 2 . General data for 10 selected sites. Station Name PGArefdownhole [m/s 2 ]

 2

			Vs30 [m/s]	Depth [m]
	HRSH13	3.08E-01	399	200
	IWTH21	5.03E-01	521	100
	IBRH16	7.05E-01	626	300
	AKTH16	1.45E+00	375	154

Table 3 . Relative comparison between the main peak frequencies in the observed BSR, the estimated BSR, and BSRlinear.

 3

		FKOH03	OITH11	KMMH03
	BSRObserved to p	2%	8%	2%
	BSRObserved to BSRlinear	15%	58%	50%

Acknowledgments

All authors are thankful with the National Research Institute for Earth Science and Disaster Resilience (NIED) in Japan for making Kik-net data available and collected the used data in this paper. Those data can be obtained from web site www.kyoshin.bosai.go.jp (last accessed October 2017).

We are thankful to the anonymous reviewers for their fruitful comments that improved the manuscript. Also, we thank the professor Fernando Lopez Caballero from Centrale Supelec, Paris-Saclay France, for his help to interpret the CyberQuake analysis.

The data in this study for any station and ground motion (e.g. fsp, PGA, PGV) are available upon request.

This research has been funded by a PhD Grant and work environment from CEREMA.

We then assume that m can be used as a borehole transfer function for the specific considered earthquake.

With this assumption, using m and the ground motion at downhole it is possible to compute the ground motion at the surface :

Where Asurface represents the Fourier transform of the accelerogram at the surface, Adownhole is the Fourier transform of the accelerogram at downhole, and m is the estimated borehole spectral ratio computed with Equation (7).

4 Results

Application to the whole KiK-net database

We compute the fsp curves for the 650 considered sites (Figure 10). In order to take into account only the sites where a reliable fsp could be derived, we choose to select the 462 sites having a standard deviation lower than 0.07. On figure 10, the curves are drawn with plain lines in the PGA range covered by the data and are represented with dashed lines in the PGA range where the curve is extrapolated. At some sites, the curve is linear over the whole PGA domain, meaning that there is no shift of frequency no matters the intensity of the ground motion is. In other sites, the decrease of fsp starts at low PGA, meaning that even for weak motions a shift of the peak frequencies is significant. The variability of the curves reflects the variability of the non-linear soil behavior at the different sites. The figure shows as well, that for some PGA at the downhole, all of the sites are prone to develop non-linear soil behavior. For example, for the evaluated stations in this study, 50% of the stations have an fsp lower than 0.8 for a PGAdownhole of 4.4x10 -1 m/s 2 . 25% have developed a lower value of fsp for a PGAdownhole of 2.3 x10 -1 m/s 2 , and 75% for PGAdownhole of 8.4x10 -1 m/s 2 . Among all of the KiK-net sites, we performed a selection of 8 sites to show in more details the fsp curves. We select sites with different range of PGArefdownhole and sites that have recorded very large earthquakes as Fukushima or Kumamoto earthquakes. The fsp curves derived at these stations are already shown in orange in Figure 10. They are presented in Figure 11 together with the individual fsp obtained for each recorded earthquake. In this figure, the graphs are organized from very non-linear stations to linear ones. The Table 2 shows the PGArefdownhole, the Vs30, and the depth where the downhole station is located.


Ibrahim Khalil
email: ibrahim.khalil@ensicaen.fr

Carlos Celis Cornejo

Karine Thomas

Philippe Bazin

Arnaud Travert

Víctor David J Pérez-Martínez

Gabriel Baldovino-Medrano

Jean-François François Paul

Françoise Maugé

Carlos M Celis-Cornejo

Víctor G Baldovino-Medrano

In situ IR-ATR study of the interaction of nitrogen heteroaromatic compounds with HY zeolites: experimental and theoretical approaches

Keywords: IR-ATR Infrared Spectroscopy -Attenuated Total Reflectance DFT Density IR-ATR in situ, organonitrogen compounds, liquid-solid interface, indole, pyridine, zeolite, Brønsted acidity, silanol group, molecular modeling

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Table of contents

The interaction modes leading to inhibit the hydrocracking catalysts (e.g. zeolites) by the weakly basic nitrogenates, such as indole, were not yet fully understood. The low vapor pressure of these molecules was presenting the inherent difficulty to their study in the gas phase setups. In this work we evidenced, by combining liquid-solid IR-ATR experiments with Density Functional Theory calculations, the vibrational frequencies of the different interaction modes of pyridine and indole with the zeolitic acid sites.

INTRODUCTION

Understanding the adsorption of nitrogen-heteroaromatic compounds over zeolites is a key to improve their design and applications in the chemical industries. [1] The interaction of nitrogen heteroaromatics with zeolites has been amply studied within the frame of refining processes such as hydroprocessing and hydrocracking. [2][3][4][5][6] Among these compounds, pyridine and indoletype molecules are commonly present in the petroleum feedstocks of the refineries. [7] They are recognized to exert inhibitory or poisoning effects over the activity and selectivity of zeolitebased catalysts employed in hydrocracking units. [8][9][10] The inhibitory effect induced by strongly basic pyridine type compounds; which possess an electron-pair donor ability, is broadly documented. [10][11][12] Nonetheless, the effect of nitrogen heteroaromatics such as indoles (pKa Indole-H + /Indole = -3.6), carbazoles, benzocarbazole over hydrocracking catalysts has only been recently put in evidence. [4,9] Indeed, a previous work confirmed that weakly basic nitrogenates are the most abundant nitrogenates in a hydrotreated vacuum gas oil (ca. 20 ppm of total nitrogen content). [4] Even at such a low concentration level, these compounds reduced the catalytic activity and influenced the selectivity of the employed catalyst. [4] However, the interaction mode leading to inhibit the zeolitic acid sites by the weakly basic nitrogenates was not fully understood.

Otherwise, N-heteroaromatic compounds such as indole are also valorized in the synthesis of fine chemical (pharmaceuticals, agrochemicals, and dyes). [13][14][15][16][17][18][19][20] These reactions are conventionally performed over both liquid phase Brønsted acids and Lewis acids. [21] However, due to the production of harmful wastes, corrosion issues, and the difficulty of product separation, recent efforts have been made in replacing these substances by zeolites. [19,20] Indeed, the reaction of indole with zeolitic acid sites was reported for the production of biologically and physiologically active chemicals, by electrophilic or nucleophilic substitutions reactions. [22] Gopal et al. studied the vapor-phase alkylation of indole with methanol to produce 3methylindole over a series of protonic and rare earth cation exchanged Y zeolites and HZSM-5 zeolites. [19] The authors reported that the Brønsted and Lewis acid sites of medium strength are key for the performed reaction. Karthik et al. investigated the electrophilic substitution over the C3 (carbon 3) atom of indole with aldehydes, using proton exchanged HHY, HZSM-5, and Zn ion-exchanged Y zeolites, for producing bis(indolyl)methanes; that are active substances for the promotion of beneficial estrogen metabolism in the human body. [20] The authors claimed a correlation between the increases of the reaction yield and the amount of Lewis acid sites, assessed by DRIFT measurements of adsorbed pyridine. However, the interaction mode of indole molecules with the active sites in the zeolites was not investigated.

In this sense, IR spectroscopy can be used to discriminate the different interaction modes of the molecules with the acid sites of the catalysts. [23] Most of the work, in which IR spectroscopy characterizations are included, make use of the transmission-adsorption and diffusive reflectance mode of analysis. [24] Thus, IR cells for in situ and operando studies working on transmission mode have been developed for characterizing the different interaction modes of the molecules with the acid sites after specific activation conditions and/or during the reaction. [23] The feasibility of an IR gas phase study is limited by the vapor pressure of the studied probe molecules. For pyridine, the saturation of solid materials under low pressure conditions is relatively easy to perform since the vapor pressure of the molecule is about 2733 Pa at 298K. [25] However, when the molecules exhibit low vapor pressures, such as indole and indolederivatives (P sat indole = 1.6 Pa at 298K), their handling in the IR in situ cell becomes either problematic or unfeasible. In such cases, the characterization with IR-ATR spectroscopy opened a wide range of possibilities for in situ studies of liquid-solid adsorption. [26] Experimental setups concerning IR-ATR continuous flow, [27][28][29] batch cells, [30] and others have been adapted for studying the adsorption of probe molecules in the liquid phase, [26] as well as for homogeneous and heterogeneous catalytic processes and electrochemical studies. [31][32][33][34][35][36][37][38] In particular, Rivera et al. performed pyridine adsorption in liquid phase over sol-gel films deposited over the ATR crystal solubilized in n-heptane and methanol. [27] Panella et al. used a liquid-solid IR-ATR batch cell to study the reactivity of amino-functionalized silica-coated magnetic nanoparticles of SiO 2 /Fe 3 O 4 . [30] As a complement to IR spectroscopy experiments, first-principle calculations for the adsorption of probe molecules can provide deeper insights into the mechanisms of adsorption. Figure SI-1 depicts the formation of the ion-pair complex for pyridine adsorbed on a model Brønsted site. [39- 41] In this sense, Castellà-Ventura et al. described the vibrational modes of pyridine adsorbed on a twelve tetrahedra (12T) cyclic cluster, for modeling the zeolite acid site. [42] The authors drew attention to the importance of including dispersive interactions in the quantum mechanical model for future studies. Their results were in line with the formation of an ion-pair complex of pyridine H-bonded to the acid site, concluding that the mentioned adsorption mode is energetically more stable than the non-protonated hydrogen-bonded pyridine molecule. Thus, the interest of the theoretical assessment for predicting the characteristic bands in the IR spectrum for a specific interaction of probe molecules was clearly shown. We should mention here that the literature does not report any study of the adsorption of indole molecules on acidic solids.

Given the foregoing, in this work, an assessment of the adsorption of pyridine and indole on zeolites via an in situ liquid-solid IR-ATR methodology is presented. Two zeolites presenting different Si/Al ratios were studied: namely, HY and USY. Besides, siliceous MCM-41 was also used as a material free of aluminum and with very limited acidic properties. Experiments were complemented by DFT calculations to help in the interpretation of the results.

In the first part of the study, the validation of the IR-ATR liquid system was successfully performed by comparing pyridine adsorption in the gas and liquid phases over different solids.

The bands relative to the different adsorption modes of pyridine, e.g. over Brønsted acid sites, silanol groups, and Lewis acid sites were distinguished from the characteristic bands belonging to pyridine in the liquid phase. In the second part, the adsorption and protonation of liquid indole, an unfeasible experiment in the IR gas setups, was successfully carried out in the IR-ATR liquid-solid setup. IR-ATR allows also exploring the low wavenumber zone (<800 cm -1 )

which is generally opaque in classical IR gas analyses, mainly due to the presence of structural vibrations of the studied materials. Combining experiments and DFT calculations provided evidence for the protonation of indole. Thus, the obtained results confirm that the weakly basic nitrogenates such as indole may strongly interact with zeolitic Brønsted acid sites, which explains the inhibitory effect on acid-type reactions reported on zeolite catalysts. [4,9] 2. EXPERIMENTAL SECTION

Studied Materials

Three solids were studied: namely, proton exchanged HY zeolite, Si/Al = 2.5, ultra-stable Y (USY) zeolite, Si/Al = 33, and a siliceous MCM-41. HY zeolite and H-USY were supplied by Union Carbide and Zeolyst International, respectively. Both zeolites were used as provided by the suppliers. MCM-41 was prepared using the procedure of Grün et al., as described in Supporting Information (SI). [43]

Assessment of Physicochemical Properties

Chemical composition. The silicon and aluminum contents of the zeolites were measured by inductively coupled plasma (ICP) coupled with optical emission spectroscopy using a Varian ICP-OES 720-ES. ca. 50 mg of the solid to analyze was first dissolved in 1 ml of aqua regia (HNO 3 / HCl = 1/3 v/v) and 3 ml of hydrofluoric acid (Merck, ≥48% in water). After 1 h of stirring at 383K, 20 ml of ultra-pure water (VWR, > 99.99%) and 3.6 g of boric acid (Merck, ≥99.5%) were added. Finally, the solution was diluted until reaching a total volume of 100 ml.

The concentrations of Si and Al obtained from ICP-OES are shown in Table 1.

Porosity and surface area. The porosity and surface area of the solids were estimated from nitrogen adsorption-desorption isotherms measured at 77K using an ASAP 2020 (Micrometrics)

gas adsorption system within a relative pressure (P/P 0 ) range between 0.05 and 1.00. Before the adsorption of nitrogen, samples of ca. 100 g were outgassed at 623K under vacuum for 4 h. The total and external surface areas were determined based on the BET theory for the zeolite samples and on the -plot method for the mesoporous silica solids, [44,45] using non-porous silica LiChrospher Si-1000 (surface area = 26.2 m 2 .g -1 ) as a reference. The corresponding isotherms are shown in supporting information (Figure SI-2).

IR study of the gas probes over the solids. The evaluation of the acidity of the solids was performed by gas-phase pyridine adsorption followed by Fourier Transform Infrared (FT-IR) spectroscopy. Samples were pressed into self-supported wafer (10 -15 mg, precisely weighted) with a surface of 2.0 cm 2 under a pressure of ~10 7 Pa using a Specac press. These wafers were in situ activated under secondary vacuum from 298K to 623K (1 K.min -1 ) and kept under the latter temperature for 4 h under 10 -4 Pa. FT-IR spectra were recorded after submitting wafers to the latter treatment. The vapor of pyridine was introduced into the cell at 298K and dosed from 0.09 to 0.90 µmol. A final equilibrium pressure of 266 Pa was established in the IR cell.

Desorption was carried under secondary vacuum (10 -4 Pa) at different steps of temperature until 423K for 20 min in each step. FT-IR spectra were recorded with a Thermo Fischer 6700 instrument equipped with an MCT detector. 64 scans were accumulated for each measurement with a resolution of 4 cm -1 . Spectra were analyzed by subtracting spectra collected after adsorption (or desorption) from the spectra corresponding to the activated samples. All spectra were normalized to a disc constant mass of 5.0 mg.cm -2 of the dried solids. The amount of acid sites was calculated by integrating the bands related to the residual pyridinium ions (Brønsted sites) and coordinated pyridine species (Lewis sites) after desorption and considering their molar absorption coefficients. [46]

IR-ATR Study of the Liquid Probes over the different Materials

The adsorption of solutions of n-hexane with pyridine (1 wt.%) or indole (0.5 wt.%) over the different zeolites and MCM-41 films was studied by in situ IR-ATR measurements.

The schematics of the home-built continuous flow IR-ATR setup are shown in Figure 1. The cell consists of a dome-like device mounted on the IR-ATR crystal that allows the flow of either gases or liquids streams by a 3-ways valve. The Internal Reflection Element (IRE) was provided with a Specac Golden Gate ATR device along with a diamond crystal and one reflection beam. The IRE was placed in a Nicolet iS50 FTIR spectrometer equipped with nitrogen cooled MCT detector and purged under a flow of dry air without CO 2 . Two recipients, one containing the solution with pyridine or indole and another containing the pure solvent were connected simultaneously to a multi-channel peristaltic pump (Ismatec). An additional 3-ways valve was employed for switching from the flow of the solution to the solvent. The solvent employed was n-hexane (99%, Sigma-Aldrich), and the test solutions contained 1 wt.% of pyridine (99.9%, Merck) or 0.5 wt.% of indole (>99%, Sigma-Aldrich). The indole solution was less concentrated due to the low solubility of this compound in non-polar solvents. Finally, a recipient was employed for storing the liquid effluents at the outlet of the IR-ATR cell. Prior to adsorption tests, a background was always recorded with the empty cell under an Argon flow of 50 ml.min -1 . Each spectrum corresponds to the average of 32 scans collected with a resolution of 4 cm -1 (time resolution 20 seconds). Spectra from the liquid solution and the solvent were recorded. The samples of the zeolites were pretreated ex situ at 623K for 4 h under secondary vacuum (10 -4 Pa). The method proposed by Davantès et al. was adapted to deposit the solid film of the sample onto the ATR crystal. [28] Accordingly, a droplet of the concentrated suspension of the pretreated solid, with a granulometry less than 20 µm, was put on top of the IR-ATR crystal and left to dry. Thereafter, the deposited zeolites were dehydrated in situ under an argon flow of 30 ml.min -1 at 473K for 4 h. The MCM-41 film was studied without in situ activation because of its instability during the treatment. After activation, the system was cooled down to RT, and an n-hexane flow was put in contact with the solid sample, we should mention that, in this step, all the solid films fairly resisted the sheer tension of the liquid flow. From this

Pump Solution Solvent

Flowmeter Ar Continuous flow ATR-IR Cell point onwards, a loop for the acquisition of IR-ATR spectra every 20 s was used to follow the accumulation of n-hexane. The spectrum corresponding to saturation of the solid film by nhexane was taken as a reference to perform further difference with the spectra corresponding to adsorption of the nitrogen aromatic compounds in solution with n-hexane on the solid film.

Afterward, keeping the same flow rate, the flow was switched to the solution containing pyridine or indole until reaching saturation. Finally, the desorption of the weakly or nonadsorbed probe molecules was achieved by switching back the liquid flow to n-hexane for several minutes.

DFT Calculation of the Frequencies of Organonitrogen Compounds on

Protonic sites

The molecular level adsorption of the nitrogen heterocycles was assessed by theoretical DFT calculations performed on a Brønsted acid site located at the pore mouth of a model mordenite. [47] Mordenite was chosen instead of faujasite to save computational resources, because of the considerably reduced number of atoms in the conventional unit cell of the latter;

144 atoms for mordenite in comparison to 576 atoms for faujasite, if we compare the pure siliceous structures. In addition, the modeling was a constraint to the adsorption of the organonitrogen compounds on the pore mouth of mordenite. In this regard, the choice of mordenite, which has a 12T main cavity as in the supercage of faujasite, should not significantly affect the calculated vibrational frequencies. Such an assumption is based on a theoretical finding that evidenced that the stretching mode of the O-H group of Brønsted acid sites is not significantly affected by the position of the aluminum atom in the different cavities of the framework. [48] Additionally, the IR spectra of pyridine adsorption on different Y-zeolites and mordenite showed the characteristic bands of pyridine almost at the same frequency values. [49] It was also shown that the strength in the acidity, which is well described by the adsorption energy of electron donor basic molecules, is more related to the chemical composition of the zeolite than to its morphology. [50] Calculations were carried out using the Vienna Ab-initio Simulation Package (VASP) considering periodic boundary conditions and pseudo potentials. [51][52][53] The Projected Augmented Wave (PAW) method and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional were chosen for performing all the calculations. [54][55][56] An energy cut-off of 600 eV was set with a -point for sampling the Brillouin-zone. Gaussian smearing with a sigma parameter of 0.05 was chosen. Grimme's density functional dispersive correction with zero dampings (DFT-D3)

was included to account for dispersive interactions. [57,58] The external surface of the (001) plane of a mordenite zeolite was modeled following the procedure described in previous work. [59] The Brønsted acid site at the pore mouth was simulated by replacing one silicon atom for an aluminum one, and subsequently hydrogenating the adjacent oxygen atom for compensating the charge. The T2 position of the Brønsted acid site on the twelve tetrahedra principal cavity was chosen following Demuth et al. assumption that this represents the most stable configuration. [60] Terminal silanol groups were created when completing with hydrogen atoms the dandling oxygen atoms bonds that are generated when the 001 plane is cut. A vacuum of 10 Å in the zaxis is added. The optimized 193-atom structure of the external surface of mordenite is shown in supplementary information (Figure SI-3). The acidic characteristics of the modeled Brønsted site in the pore mouth of the mordenite were validated in previous work, [59] calculating the vibrational frequency of the O-H group adjacent to the aluminum atom, and the atomic charge of the hydrogen atom. A charge of +0.65 and a frequency of 3604 ± 5 cm -1 was taken to represent the characteristics of a Brønsted acid site. [48,61] For studying the adsorption of the probe molecules, pyridine and indole were first optimized in a vacuum cube box with an edge length of 25 Å and a frequency calculation was subsequently performed for obtaining the data of the isolated gas-phase molecules. The predicted vibrational frequencies (not scaled) for pyridine and indole in the gas phase, interacting with the silanol groups or with protonated OH groups of the zeolite were calculated. The interaction of pyridine and indole over a Lewis acid site was not computed in this work. The protonation of indole was performed over the carbon in the position 3 of the molecule (C3 site), known and checked to be the site with the highest proton affinity (PA) from all the 7 sites where the protonation is possible (PA (C3 protonation) = 212.4 kcal/mol). [62] 3. RESULTS AND DISCUSSION

Chemical and Textural Properties of the Solids

Table 1 presents the chemical composition and textural properties of the different materials.

From ICP results, the Si/Al ratio of the zeolites were calculated. An increase in the Si/Al ratio for the zeolite USY33 containing less aluminum than HY2.5 is observed. The mesoporous volume of the USY zeolite is higher in comparison to the HY2.5 zeolite due to the partial destruction of the structure during the aluminum extraction procedure performed by the supplier. The measured chemical and textural properties of the zeolites are similar to the characterization results obtained for these materials in the literature. [63][64][65][66] The MCM-41 solid did not present any microporosity as shown in Table 1 and Figure SI-2 (nitrogen sorption isotherms), also the obtained surface area, calculated using the -plot method, was around 1000 m 2 .g -1 which is in a good agreement with the results shown by Grün et al. [43] Table 1. Chemical composition and textural properties of the different materials

Adsorbents

ICP results (mg.g -1 ) Surface area a (m 2 .g -1 ) Pore volume (cm [45] 3.2 Pyridine Adsorption

IR Study of the Adsorption of Gas phase Pyridine over the Solids

In a first step, the acid sites were studied by pyridine gas phase adsorption on the solids activated in situ at 623K. The IR spectra of the adsorbents before pyridine adsorption are shown in Figure SI-4 and commented in the Supplementary Information. Figure 2 presents the IR spectra of the pyridine species remaining after in situ desorption at 423K under secondary vacuum. The results show the presence of two interaction modes of pyridine with the acid sites.

The band at 1545 cm -1 (v 8a mode) is characteristic for pyridinium ions PyH + formed over Brønsted acid sites, while the band at 1455 cm -1 (v 19b mode) reveals coordinated pyridine (LPy)

over Lewis acid sites (Figure 2). [63,67] The observed band at 1490 cm -1 is characteristic of the v 19a mode over both Brønsted Acid Site (BAS) and Lewis Acid Site (LAS), although its sensitivity to the concentration of BAS is higher. [68] The concentration of both BAS and LAS can be calculated from the integration of v 8a and v 8a bands, respectively, using the following absorption coefficients:  PyH  ) = 1.8 cm.µmol -1 and  (LPy) = 1.5 cm.µmol -1 . [46] Figure 2. IR spectra of the residual pyridine species over the different solids after desorption at 423K under secondary vacuum. 2 shows that the amount of BAS strongly decreases on USY33 as compared to HY2.5, whereas, the amount of LAS increases due to the formation of extraframework aluminic phase species during the dealumination treatment. [69,70] Over MCM-41, which only presents weakacidic silanol sites, no acid sites were evidenced after the pyridine desorption at 423K. These results confirm that the studied solids present various properties in terms of nature and concentration of acid sites that can serve as references to validate the characterizations developed herein.

IR-ATR Study of Liquid phase Pyridine Adsorption over the Solids

IR-ATR study

IR-ATR spectra of the thin films of Y zeolites and MCM-41 after deposition are shown in Figure 3. It must be underlined that for IR analysis in transmission and diffusion modes, of a self-supported silica-based wafer, the structural vibration bands of the solids make the zone below 1200 cm -1 almost completely opaque. [69] By contrast, as shown in Figure 3, these bands present limited intensities in the IR-ATR spectra of thin films, which allows their analysis. For the zeolite samples, the bands in the region at 850-650 cm -1 are characteristic of the mix of the symmetric stretching and bending modes of the T-O-T bridges (2). Also, these bands shift up, the lower is the Al content, as shown for the bands at 789 and 724 cm -1 over the HY2.5 zeolite that shift up to 831 and 787 cm -1 over the low aluminum loaded USY33 zeolite. [69] However, it should be noted that probe molecule adsorption can induce distortion at the surface of the solids that can shift the wavenumbers of these structural bands. Thus, positive -negative bands can appear on the difference spectra. Hence, caution must be taken for the interpretation of the bands detected in this low wavenumber zone after the adsorption of pyridine or indole.

In general, in IR-ATR mode, all the IR bands as well as those corresponding to the adsorbed species present very low intensities due to the limited thickness of the samples that are analyzed (a few microns). Thus, cautions must be taken in order to avoid any artifacts that can hide the bands coming from adsorption. Among them, one may mention that it is important to watch the stability of the concentration of water present as a contaminant in the purge of the spectrometer during the experiment. Special care must be taken for spectral analysis of the region between 1400 and 1600 cm -1 where IR bands of water are very intense.

IR-ATR study of the liquid phase pyridine -solid interaction

IR-ATR spectra of the liquid pyridine adsorption over the films of HY2.

High wavenumber zone analysis (1650-1400 cm -1 )

The spectrum of adsorbed pyridine on the HY2.5 zeolite (Figure 4A -solid line) presents bands at 1544 and 1490 cm -1 , which characterize pyridinium ions formed over BAS. These bands remain after flushing with the solvent (Figure 4A -dashed line), hence confirming that they correspond to strongly adsorbed species. The bands at 1595 and 1444 cm -1 were attributed to weakly H-bonded pyridine molecules (Py---H) with the solid, [71,72] this adsorption mode did not resist the desorption under n-hexane flow. The bands for liquid pyridine were also detected during the adsorption (1578 and 1437(shoulder) cm -1 ), however, they were also eliminated during the desorption.

Adsorption results obtained over USY33 are shown in Figure 4B. The bands related to pyridinium ions (Figure 4B -solid line) are present at similar positions than over the HY2.5 zeolite (1544 and 1490 cm -1 ). Note that the intensity of the BAS bands is less than on HY2.5 as expected from the lower Al amount of USY33 (Table 1). The bands at 1597 and 1446 cm -1 indicate the presence of LAS on this partially dealuminated zeolite (Table 2). [63] This assignment is confirmed by the stability of these species under flowing hexane. In addition, the bands of pyridine interacting with silanol groups (1596 and 1444 cm -1 ) and liquid phase pyridine (1578 and 1437 cm -1 ) were also detected. However, these bands disappeared by flushing with n-hexane (Figure 4B -dashed line).

Over the MCM-41 film, pyridine bands are not intense (Figure 4C -solid line), even being hard to tell them apart from the spectral noise. The spectra mainly exhibit the bands for liquid pyridine. Only two bands at 1596 and 1444 cm -1 indicate the formation of weakly-adsorbed Hbonded pyridine species. Rivera et al. reported a similar pattern for pyridine adsorption in the liquid phase on silica, [27] which allowed observing mainly the characteristic bands of liquid pyridine. They reported the appearance of one weak band at 1492 cm -1 attributed to pyridinium ions, indicating the presence of strong acid sites. This suggests some Al contamination in the silica. These authors also highlighted the presence of noise in the spectra ascribed to water in the vapor phase that distorts the spectra between 1600 and 1400 cm -1 , which makes difficult to distinguish pyridine bands in this region.

The acidic site nature of the solids detected by pyridine under liquid phase is in good agreement with that obtained with pyridine under gas phase (section 3.2.1 -Table 2). It is not possible to make a quantitative assessment of the acidic sites by IR-ATR due to the absence of accurate knowledge of the volume of solids probed by the ATR signal. However, a comparison of the ratio between the various pyridine species formed over the tested solids in the gas and liquid phases experiments can be done. The ratios of the areas of the bands located at 1490 cm -1 (Hbonded and coordinated pyridine) and 1544 cm -1 (protonated pyridine) were compared for pyridine adsorption performed in the gas-solid conditions in the IR transmission set up (areas measured after pyridine evacuation at 423K) (Figure 2) and performed in the liquid-solid with IR-ATR set up (areas measured under hexane flow) (Figure 4). Note that the widespread of the noise across the IR-ATR spectra limits an accurate measurement of the band area. The ratios reach values of 1.3 (gas experiment) and 1.5 (liquid experiment) over USY33. Over HY2.5 zeolite, the ratios for gas and liquid experiments were, even more, close with a value of almost 1.1 in both experiments. A good agreement between gas and liquid phase adsorption of pyridine is thus shown.

Low wavenumber zone analysis (800-650 cm -1 )

The study of the zone from 650 to 800 cm -1 presents the great advantage to overcome the problem of the spectral noise induced by the water present in the vapor purge. Indeed, water vapor does not present any vibration-rotation bands in this zone and IR-ATR allows the analysis of this spectral zone for thin films of zeolites and oxides. However, as mentioned previously, some solid structural bands are detected in the region of low-frequency (800 -650 cm -1 ) (Figure 3) and these bands can be shifted by structural distortion due to pyridine adsorption on the solid.

This warns that caution should be taken for the interpretation of the bands in this zone. As for the HY2.5 zeolite, an assignment of bands characteristic of weakly adsorbed species or nonadsorbed ones can be proposed by comparing pyridine spectra obtained before and after flowing hexane regarding the spectrum of liquid pyridine (Figure 4A -right). Hence, the band at 706 cm - 1 can be attributed to H-bonded pyridine, while the bands at 743 and 703 cm -1 are related to pyridine in the liquid phase. The strong band at 686 cm -1 that resists to hexane flushing likely accounts for pyridinium species. However, due to the presence of strong (str) structural bands between 800-740 cm -1 on HY2.5 (Figure 3), the bands detected at 787 (str), 760 (str), 728 (wweak) cm -1 after pyridine adsorption can arise from a shift of the structural band or pyridine species. Their assignment being shaky with only experimental data, pyridine adsorption frequencies were computed by DFT calculation to validate the proposed assignments (see next section).

On USY33 zeolite, the structural bands present limited intensity in the region below 780 cm -1 (Figure 3), allowing an easier assignment of the pyridine adsorption bands on this solid. During the adsorption of pyridine, new bands appear at 750(str), 727(w), 703(str), 686(m-medium) cm -1

(Figure 4B and its zoom in Figure 5). The intense bands at 750 and 703 cm -1 , that are partially eliminated by flushing, can be assigned to pyridine adsorbed over silanol groups. This is in good agreement with the assignment proposed in the previous section. The bands detected at 727 and at 686 cm -1 whose stability under flushing is clearly observed in Figure 5, should characterize protonated pyridine as proposed on HY2.5 zeolite. A more robust assignment will be done after comparison with the frequencies calculated by DFT (section 3.2.3).

In summary, a good agreement was observed between pyridine adsorption in gas and liquid phases. However, IR-ATR spectra present low intense bands, and the analysis of the region between 1600 and 1400 cm -1 was difficult due to the presence of water bands as a contaminant in the spectrometer gas purge. [27,30] Nevertheless, an advantage of the IR-ATR setup is in giving access to analyze the low wavenumber zone (800 -650 cm -1 ), where no water bands appear.

Pyridine exhibits some intense bands in this zone that provide interesting information about the adsorption mode of pyridine in the liquid-solid IR-ATR analysis. [73] To go further, DFT calculations were performed in order to predict expected frequencies according to the interaction mode of pyridine with the surface.

Theoretical Assessment of the Adsorption of Pyridine

Most of the IR studies using pyridine for probing the acid sites focus on the region between 1400 and 1700 cm -1 . The corresponding in-plane atomic displacements for the pyridine vibrational modes appearing in this wavenumber zone ( 8a ,  8b , 19a and 19b ) are shown in Figure 6. [74] Our work points out the interest of accurately analyzing the vibrational zone between 800 and 600 cm -1 for helping in understanding the interaction between liquid pyridine and the solid films. Within the low-frequency zone,  10b and 11 modes are two characteristic signals of pyridine that correspond to out-of-plane motions (Figure 6). Table 3 gathers the frequencies in the high-wavenumber and low-wavenumber zones calculated by DFT for the various vibration modes of gas-phase pyridine, H-bonded, and protonated pyridine. The frequencies of coordinated pyridine molecules over a Lewis acid site were not computed in this work. Table 3 also compares the theoretical and experimental frequencies obtained in this work to theoretical frequencies found in the literature. [42] In general, a good agreement appears between the frequencies calculated for pyridine in the gas phase and those for pyridine adsorbed over the various sites. The calculations show that H-bonding leads to limited changes in wavenumbers compared to pure pyridine. This agrees with the fact that the interaction of pyridine with a silanol can be described as physical adsorption.

In the high-wavenumber zone, the vibrational frequencies predicted for pyridine on silanol groups are in good agreement with the values obtained experimentally (Table 3). For this kind of interaction, the characteristic vibrations of the pyridine molecules are very similar to those observed for pyridine in the gas phase. Only when the interaction is strong enough, a chemisorbed complex of pyridine on a Brønsted site, some changes are observed in the motions of the pyridine molecules. [42,75] These changes could explain the registered shifts which are characteristic of the ion-pair H-bonded complexes described in the literature. [42,76] For the interaction of pyridine with the Brønsted acid sites, the predicted shifts of the vibrational modes  8a and  19a were similar to what was previously reported in the literature. [42] However, for the  8b and  19b , vibration modes, the previously calculated values were slightly higher as compared to the experimental data. [42] The computed values obtained herein are more realistic (Table 3). Previous calculations considered a 12T OH size cluster, [42] while in the present work the adsorption occurs on a Brønsted acid site located in the pore mouth of a 193-atom mordenite.

Table 3. Comparison between the predicted vibrational frequencies and the experimental values obtained in IR-ATR for pyridine adsorption over silanol groups and Brønsted acid sites.

Pyridine frequencies (cm

Vibrational mode*

Gas phase Silanol group Brønsted site

Other † (Calc.) [42] This work (Calc.)

Other † (Calc.) [42] This work (Calc.)

This work (IR-ATR)

Other † (Calc.) [42] This work (Calc.) calculated for the pyridine interaction on a Brønsted acid site in a 12T OH zeolite cluster using the B3LYP/6-31G(d,p), as reported in reference. [42] This work presents, as far as we know, the first calculation of the vibrational frequencies of pyridine in the region of low-frequency (800 -650 cm -1 ). Two intense bands were observed:

 10b and  11 , which correspond to out-of-plane vibrations (Figure 6). For the  10b and  11 modes over a silanol, a small shift was found with bands predicted at 745 and 706 cm -1 . This agrees with the observation done for liquid pyridine adsorption over the MCM-41 film (H-bonded bands at 750 and 703 cm -1 ). Meanwhile, for the interactions with a Brønsted site, characteristic bands at 728 and 680 cm -1 of pyridinium species are predicted. [73] Hence, based on the calculated frequencies of pyridine chemisorbed on Brønsted sites, the bands observed at 727 and 686 cm -1 during the IR-ATR analysis of liquid pyridinezeolite interaction (Figure 4) can be assigned to pyridinium species.

In summary, a good agreement is observed between the predicted (by DFT calculation) and experimental (using the liquid phase IR-ATR study) vibrational frequencies of the pyridine adsorption over silanol groups and Brønsted acid sites. In the high wavenumber zone (1650 -1400 cm -1 ), as expected, the bands observed at 1544 and 1490 cm -1 were assigned to protonated pyridine molecules over the Brønsted sites of the zeolites. Most notably, in the low frequency zone (850 -600 cm -1 ) the experiments and DFT calculations allows establishing that the bands observed at 750 and 703 cm -1 correspond to H-bonded pyridine species whereas the bands at 727 and 686 cm -1 can be assigned to pyridinium ions formed over the Brønsted acid sites of the zeolites.

In the next section, the study of the adsorption and protonation of indole will be investigated over the previously studied adsorbents. To go further, DFT calculations were performed in order to predict expected frequencies according to the mode of interaction of indole with the surface.

Indole Adsorption

IR-ATR Study of Liquid phase Indole Adsorption over the Solids

liquid indole (dotted line) shows bands at 1455, 1414 cm -1 and a shoulder at 1468 cm -1 in the high wavenumber zone as well as bands at 768, 763, 738 and 716 cm -1 in the low wavenumber zone. The solution flowed over the ATR crystal contained indole diluted into hexane (0.5 wt.%). To obtain the indole spectrum, the bands of hexane must be eliminated by subtracting the spectra of pure n-hexane. Since n-hexane presents intense bands at 1468 and 728 cm -1 , the shoulders observed in the indole spectrum at these wavenumbers should correspond to partially eliminated n-hexane bands after the subtraction. Thus, these two bands will not be considered in the assignment of the indole bands.

After indole adsorption over the HY2.5 zeolite film (Figure 7A), bands appear at 1617, 1608, 1543, 1489, and 1455 cm -1 , in the high-wavenumber zone. In the low-wavenumber zone, bands are detected at ~758, 738, 716 and 702 cm -1 . However, a broad massif also appears at 787-750 cm -1 . As previously noticed for pyridine adsorption on HY2. The indole adsorption over the USY33 zeolite film (Figure 7B) shows new bands at 1543, 1514, 1489, 1455, and 1414 cm -1 in the high-wavenumber zone and 767, 743, 718, and 702 cm -1 in the low-wavenumber zone. Most of these bands were observed on HY2.5 spectra even if with different intensities. As previously shown over HY2.5, during the desorption most of these bands lose intensity but withstand rinsing under the solvent flow. Only the band at 1414 cm -1 disappears totally which allows to attribute it to weakly-adsorbed indole species (Figure 7B and

Theoretical Assessment of the Adsorption of Indole

To assign the bands recorded on the IR-ATR spectra, DFT calculations of the vibrational frequencies of indole either in the gas phase or adsorbed on silanol or Brønsted acid sites were performed in the high and low-wavenumber zones. Indole can interact with the silanol group and with the zeolitic OH groups as depicted in Figure 9. The H-bonded interaction involves an interaction between the NH group of the indole molecule and the oxygen atom from the zeolite framework, as well as an interaction between the SiOH group and carbon of the C5 cycle. In such a configuration, the distance between the H of the SiOH group and the -carbon is about 1.95 Å. The most stable configuration of protonated indole was obtained over the carbon in position 3 (C3 site). A similar configuration was obtained by Somers et al. [62] 0.0008 a.u. The predicted vibrational frequencies for indole in the gas phase, interacting with the silanol group or with a protonated OH group of the zeolite were calculated and are reported in Table 4. For indole in the gas phase, the calculations predicted several bands in the region of high frequencies; namely, 1613, 1575, 1505, 1485, 1444, and 1416 cm -1 . These bands are attributed to in-plane vibrational modes of the 6 and 5-membered rings (Figure 10). In the region from 800 to 650 cm -1 , calculated vibrational modes corresponding to out-of-plane motions are predicted at 757, 755, 731, and 709 cm -1 .

Gas phase

Terminal silanol group Brønsted acid site Indole keeps the same atomic motions as in the gas phase and small shifts are predicted compared to gas-phase frequencies when the molecule is H-bonded to silanol groups. Table 4 shows that H-bonded indole should give rise to bands at 1615, 1573, 1496, 1482, 1442, and 1411 cm -1 in the high-wavenumber zone. Concerning the low-wavenumber zone, bands at 774, 758, 745, and 737 cm -1 should be observed. On the contrary, for indole interacting with a Brønsted acid site, the calculations predicted shifts in the C-C stretching modes of the C5 and C6 cycles. This could be attributed to the change in the vibrational modes of the adsorbed molecule, even some additional modes can be predicted over the protonated atom (C3). Hence, the supplementary band at 1554 cm -1 corresponds to a C-C stretching vibration of the protonated carbon (sp 3 ) of the indole rings (Table 4). In the region of low-frequency, the theoretical calculation of indole on a Brønsted acid site has shown a new band at 700 cm -1 attributed to outof-plane motions of the hydrogen atoms of the C5 cycle.

The vibrational frequencies predicted for gas-phase indole adsorbed on a silanol terminal group and a Brønsted site of a mordenite can be summarized as follows: (i) for the interactions of indole with silanol groups small shifts were obtained; (ii) for protonation of indole more important shifts are expected and new bands should appear. Next, the calculated vibration frequencies were compared to the experimental spectra in order to assign the observed bands and identify the indole adsorption modes on the different solids.

Comparison between calculated and observed vibrational frequencies of indole

The reference gas phase spectrum of indole (Figure SI-5, spectra from the Omnic database and Table 4) presents bands at very similar wavenumbers to the computed ones. On the experimentally IR-ATR spectrum of indole in the liquid phase (Figure 7), the most intense bands are observed at 1455 cm -1 and 1414 cm -1 , while those at 1613, 1575, 1505, and 1485 cm -1

were not detected (Table 4). This is likely due to the low amount of indole in the n-hexane mixture (0.5wt. %) and the low intensity of IR-ATR spectra.

A parallel between the computed vibrational frequencies (Table 4) and the recorded spectra of indole on HY2.5 (Figure 7A) allows assigning the bands at 1617, 1608, and 1543 cm -1 to protonated indole species. These protonated species were also detected on USY33 (Figure 7B), however, the bands were less intense which agrees with the lower concentration of BAS over the dealuminated USY zeolite (Table 2). The noisy spectra and the lower concentration of protonated indole prevent the detection of the bands at 1617 and 1608 cm -1 on USY33. As expected, the characteristic bands for protonated indole species were not detected on MCM-41.

The band at 1489 cm -1 detected on both HY2.5 and USY33 zeolites can be assigned straightforwardly to H-bonded indole (Table 4), while the assignment of the bands at 1455 and 1468 cm -1 was not as obvious. Table 4 shows that the band at 1455 cm -1 can be characteristic for both protonated indole and liquid phase indole, whereas the band at 1468 cm -1 can be assigned to protonated indole or, alternatively, to the presence of residual n-hexane in the difference spectrum, thus making difficult to assign these bands to a unique species.

In the low wavenumber zone, liquid indole presents four bands at 768, 763, 736, 716 cm -1 , which matches satisfactorily with the calculated frequencies (Table 4). The main question to answer when analyzing the low-wavenumber region is: does this spectral region confirm the formation of protonated indole species on faujasites samples? Table 4 shows that a characteristic band of protonated indole species should appear at about 700 cm -1 . First, it must and theoretical calculations. Using IR-ATR system offers two main advantages, it allows : (1) the analysis of the low frequency region (< 800 cm -1 ) that presents characteristic vibrational modes of adsorbed species, (2) studying the adsorption of low vapor pressure molecules (such as indole) that presents an inherent difficulty to be operated in the gas phase IR cells.

A parallel between the vibrational frequencies predicted by DFT calculations, and the experimental results obtained using the IR-ATR cell for studying of the solid films under continuous liquid flow shows that:

-As for pyridine, the low wavenumber zone, a newly explored zone, can also account for its adsorption modes. Indeed, the IR bands at 686 and 727 cm -1 characterize the pyridinium species formed on zeolitic OH groups, while the signals at 703 and 750 cm -1 are associated with pyridine in interaction with silanol groups.

-As for indole, the IR-ATR study points out the presence of new bands after its adsorption over Y zeolites and MCM-41. Over the zeolites, the detection of bands at 1617, 1608, 1543 and 705 cm -1 provides evidence for the formation of protonated indole species over acidic OH groups.

These findings support the fact that indole, a weakly basic nitrogenates, can be protonated over the zeolitic acid sites and thus decreasing their activity during the hydrocracking process.

Figure 1 .

 1 Figure 1. Flowsheet type scheme of the continuous flow IR-ATR cell setup. The solution recipient contains the diluted N-molecule into n-hexane; the solvent recipient contains pure n-hexane.

Figure 3 .

 3 Figure 3. IR-ATR spectra of the thin films of zeolites and MCM-41 as-deposited on the ATR crystal.

Figure 4 .

 4 Figure 4. IR-ATR spectra of pyridine adsorption over films of (A) HY2.5, (B) USY33 and (C) MCM-41 in the high wavenumber zone (left hand side) and low wavenumber zone (right hand side). The solid lines correspond to the spectra taken at the steady-state adsorption of the pyridine solution; the dashed lines correspond to the spectra obtained after desorption of the weakly-adsorbed pyridine species by passing a flow of hexane at 298K and the dotted lines refer to liquid pyridine spectrum.

Figure 5 .

 5 Figure 5. IR-ATR spectra of pyridine adsorbed on USY33 and comparison with pure pyridine (zoom of the Figure 4B right). The solid line corresponds to the steady-state adsorption of the pyridine solution; the dashed line corresponds to the spectrum obtained after desorption of the weakly adsorbed species by passing a flow of hexane at 298K and the dotted line refers to liquid pyridine spectrum.

Figure 6 .

 6 Figure 6. Atomic displacements for the commonly studied vibrational modes of pyridine.

5 (

 5 Figure 4), this massif was attributed to the change in the frequency of the vibrations of the T-O-T bridges during the indole adsorption. During desorption under n-hexane flow, most of these bands resist flushing and remain present with lower intensities (dashed line).

Figure 8 )Figure 7 .

 87 Figure 8). The zoom presented Figure 8 makes clearly visible that the bands at 728, 718 and 706 cm -1 resist to the n-hexane flush.

Figure 8 .

 8 Figure 8. IR-ATR spectra of the indole adsorption over USY33 zeolite and of pure indole (Zoom of the Figure 7B -right). The solid line corresponds to the spectrum of the steady-state adsorption of the indole solution; the dashed line corresponds to the spectrum obtained after desorption of the weakly adsorbed species by passing through the solvent at 298K and the dotted line refers to indole spectrum in n-hexane. [* this component can correspond to a characteristic band for n-hexane]

Figure 9 .

 9 Figure 9. The three evaluated configurations for the calculation of the vibrational frequencies for indole adsorption.

Figure 10 .

 10 Figure 10. Atomic displacements for in-plane vibrational modes of indole in the gas phase in the region of 1400 to 1600 cm -1 .

Table 2 . Concentration of acidic sites in the different materials

 2

	Adsorbents	Acid sites (µmol.g -1 ) BAS LAS Total
	HY2.5	1263	22	1285
	USY33	98	61	159
	MCM-41	0	0	0

  Wilson's numbers for the vibrational modes of the benzene ring. † Theoretical data from other work,

									This work
									(IR-ATR)
	8a	1595	1587	1612	1595	1597	1634	1634	-
	8b	1591	1574	1591	1573	1578	1669	1621	-
	19b	1445	1432	1464	1437	1444	1563	1547	1544
	19a	1486	1472	1487	1475	1485	1486	1485	1490
	11	-	742	-	745	750	-	721	728
	10b	-	699	-	706	703	-	686	686

*

Table 4 . Comparison between the predicted vibrational frequencies and the experimental values obtained in IR-ATR of the indole adsorption over a silanol group and a Brønsted acid site.

 4

	Vibrational			Indole frequencies (cm	-1 )
	modes	Gas phase	Liquid phase	Silanol group		Brønsted acid site	b
		[Ref] a	This work (Calc.)	This work (IR-ATR)	This work (Calc.)	This work (IR-ATR)	This work (Calc.)	This work (IR-ATR)
		1614	1613	-	1615	1617	1623,1600 1617,1608
		1576	1575	-	1573	-	-	-
	In-Plane Quadrant and	-	-	-	-	-	1554	1543
	Semicircular	1505	1505	-	1496	1489	-	-
	Stretching of the rings	1486	1485	-	1482	1489	1468	1468
		1456	1444	1455	1442	1445	1450	1455
		1415	1416	1415	1411	1414	-	-
		765	757	768	774	767	768	768
	Out-of-plane	747	755	763	758	-	740	738/746
	vibrations of the rings and	731	731	738	745	743	717	716
	H atoms	719	709	716	737	-	-	-
		-	-	-	-	-	700	705

a OMNIC database, Aldrich Catalog No: I-340-8 b C3 protonated indole

ACKNOWLEDGMENTS

The authors thank Yoann Levaque for helping in assembling the IR-ATR setup. Labex EMC3, ANR and FEDER are also acknowledged for the PhD grant of IK and the financial support to the BIOCAR project. Authors thank Ecopetrol S.A. for funding the present work through the project AC2-5211770. The gentle collaboration of the SC3 work team from the GUANE platform at UIS is acknowledged. VGBM thanks VIE-UIS and Colciencias are also thanks for funding and for PhD grant of CMCC. The authors also thank Lille University Calculation center (CRI) partially funded by FEDER for CPU allocation.

be underlined that the experimental spectrum of liquid indole or indole on MCM-41 does not present any band at about 700 cm -1 (Figure 7C). On HY2.5, and similarly on USY33, Figure 7A and Figure 7B (and more clearly the zoom presented in Figure 8) point out a weak band at 702 cm -1 which shifts to 705 cm -1 (706 cm -1 for USY33) after flowing with the solvent. Thus, both high and low-wavenumber zones provide spectroscopic evidence for indole protonation through the appearance of bands at 1543 cm -1 (as well as 1617, 1608 cm -1 ) and 705 cm -1 . The other features detected at 767, 746-738, and 716 cm -1 both on HY2.5 and USY33 can support indole protonation but can be also assigned to H-bonded indole.

Over the MCM-41 film (Figure 7-C), only the liquid indole bands were detected in the high wavenumber zone. In the low-wavenumber zone, a small shoulder was observed at 745 cm -1 attributed to weakly adsorbed indole molecules over a silanol group. All bands disappear after flowing with n-hexane, which confirms the weak interaction of indole with the siliceous film.

Considering the collected evidence, this IR-ATR study reveals that indole, a weakly basic compound, can be protonated on acidic Y zeolites as unambiguously evidenced by the formation of the bands at 1617, 1608, 1543, and 705 cm -1 . In this way, Laredo et al. attempted to separate basic from weakly basic nitrogenates in a heavy oil fraction by using solid acids. [77] The separation was not efficient, and this was ascribed to the strong adsorption of the weakly basic compounds on the employed materials. Therefore, the interaction of indole with Brønsted acid sites via the formation of an ion-pair complex is feasible, contrary to the suggested mechanism of dissociative adsorption. [78] Such interactions should also alter the surface chemistry of these materials and thus block the Brønsted acid sites. These findings draw attention to the possible role of weakly basic nitrogen compounds that can modify the catalytic performance of zeolitic materials employed for hydrocracking units. [4]

CONCLUSIONS

The interaction modes of liquid phase pyridine and indole over HY and USY zeolites and MCM-41 solid films were successfully studied by using an in situ liquid-solid IR-ATR setup