Aline Msaed

Olivier Chadebec

Jérôme Delamare

A Simple

Keywords: Magnetic Sensor, Magnetic Gradient Measurements, Identification of Electrical Currents, Reducing of Disturbance Magnetic Field

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

With the increase of the demand for power monitoring, current sensors must be more and more compact, accurate, low-cost and with wide measurement ranges. Furthermore, various principles of current measurements have been studied. The measurement of the current with a shunt introduces many problems such as Joule losses, noise and non-isolated measurements. Current Transformers (CT) have been widely used, 1 but a small (CT) has a problems with direct current measurements and a problems of consumption due to the feedback. The current sensor based on the Hall Effect presents a large linearity error. Other principles were also proposed in Refs. [2, 3], they are based on the detection of the vibration of the cable carrying the current to be identified or by the measurement of its temperature. In this paper, we are interested in the measurement of the current from point measurements of the magnetic fields.

Nevertheless, there are several difficulties with this method; Industrial environments are often situated within disturbing magnetic fields which can induce large measurement inaccuracies. At present, in a current sensor based on the measurement of magnetic fields, the major problem is to make sure that the measured field derives exclusively from the current to be identified. A single magnetic sensor is insufficient to reconstruct the current flowing in a conductor, due to the presence of the earth's magnetic field. Thus, the main problem resides in the number of magnetic sensors and the manner of their positioning around the conductors. In addition, system sensibilities to external magnetic perturbations lead to an increase in reconstructing errors especially when we have a conductor with a high current value nearby the current to be reconstructed. One solution is to shield the probes with ferromagnetic material. However, the complication size and cost of the sensors increases, 4 particularly when using double shields. Finally, the limitations of the magnetic sensors affect the measurements, such as saturation of the sensors and misalignment when using directional magnetic sensors.

We will show the results of the theoretical study for the approach of field and gradient measurements for a simple conductor by taking into count the advantages and the disadvantages for each approach. The approach of magnetic gradient measurements greatly reduces magnetic disturbance. It consists in placing non-contact magnetic field sensors nearby the conductors and deducing the currents by inverting the Biot-Savart law. 5 We will apply theses approaches for a three-phase system by illustrating the numerical results for the reconstruction of a small current in the presence of a homogeneous and non homogeneous magnetic field. But, when implementing the practical application for the gradient measurement approach, many difficulties will be appearing; we will list it in the following sections. Then, we will suppose a configuration solution able to subdue these difficulties; it is based on "roundtrip" conductors and consists to laying the two magnetic sensors collinearly along the same axis above the conductor sections to apply gradient measurements. Next, we will apply this configuration for a three-phase system by describing the matrix form and the manner in which to calculate the reconstruction of the currents.

BASIC THEORY

Let I c be a current flowing in an infinitely long conductor, a scalar magnetic sensor placed closed to the conductor measures the main tangential component of the magnetic field. The directional field sensor is placed at a distance of a few millimeters to the conductor I c ; it is represented by an arrow in (Fig. 1). Thanks to the 2D Biot-Savart law, we can easily compute the magnetic field due to I c :

H = I c 2 r = k • I c ( 1
)
Where r is the distance between the conductor and the field sensor. By inverting the previous relation, it is possible to get the current value from the measurement of magnetic field (H meas ):

I c = k -1 • H meas (2)
However, in real applications, it is usual that the magnetic sensor will be exposed to external fields of unknown origin such as a nearby conductor carrying a high current.

Let us consider an infinitely long disruptive conductor in which flows a high current I d (I d = 300 • I c ) generating an additional field: We consider that the influence of this disruptive conductor is negligible if the ratio between the fields created by I d and I c is inferior to 2%.

H meas = H I c + H I d ( 3
)
= H I d H I c < 2% (4)
Then we can plot an exclusive area where the disruptive conductor should not be (Fig. 1(a)). This figure illustrates that I d must be placed at a distance higher than 16 m around the conductor I c ; an enormous distance. This condition is impractical to enforce in real industrial applications.

The approach of gradient measurements is based on the derivative of the magnetic field. In 2D, the field generated by the conductor I c decrease as 1/r 2 for the gradient; compared to 1/r dependence for field measurements. Furthermore, the gradient measurement naturally rejects constant external magnetic perturbations.

dH dr = -I c 2 r 2 = M • I c (5)
In practice, this approach can be achieved by using two magnetic scalar sensors placed closed to each others (Fig. 1(b)) and by subtracting their two measured fields.

In this computation, the gradient measurement is considered as being punctual. If we consider the previous case, where the influence of the disruptive conductor is negligible for < 0 02, the exclusion area for disturbance conductor therefore decreases to 13 cm. Notice that this gradient-based approach can be extended to the measurements of gradients of gradients by applying the second spatial derivative to the magnetic field:

d 2 H dr 2 = 2I c 2 r 3 = L • I c (6)
Then, the field generated by I c has a 1/r 3 decrease law for gradient of gradients. It can be put into practice by using three field detectors placed close to each other as shown (Fig. 1(c)). For the above case, the exclusion area radius decreases to 3 cm (Fig. 1(c)).

PHASES GENERALIZATION

Considering three-phases, we have applied these three approaches for reconstructing a small current I c by supporting to the Biot-Savart law; n magnetic sensors surround the system according to each approach. Thanks to the linearity of Biot-Savart law, each sensor will measure the sum of the magnetic field generated by the three currents.

H 1 = k 11 • I 1 + k 12 • I 2 + k 13 • I 3 H 2 = k 21 • I 1 + k 22 • I 2 + k 23 • I 3 H n = k n1 • I 1 + k n2 • I 2 + k n3 • I 3 (7)
Then, we can write matrix system:

H i = k ij • I j i = 1 n j = 1 3 (8)
The coefficients of this matrix depend on the geometry of the conductors and the number of the field sensors and the locations. Once the geometry of the system is fixed, we can compute the matrix elements. This is done simply by imposing a known current through each conductor individually and computing the magnitude of the magnetic field from each field sensor. [START_REF] Olson | Industry Application Conferences[END_REF] This method can be applied numerically if the geometry is more complex (by using Finite Element approach, for example), by in any case, the problem should be linear. Having determined the matrix coefficient the inverse problem is solved by the least squares approach. The system to solve becomes:

K T K • I = K T • H (9)
Extracting the currents is accomplished by solving it. The approach of field measurements is applied by using one magnetic field sensor for each conductor, whereas the gradient measurements is made by using two conductors placed closed to each other upon each conductor and by subtracting their measured fields. Then, to apply the approach of gradient of gradient measurements for threephase sensors, we need three magnetic field sensors placed nearby the conductors, and the approach is executed by doing a double difference of the matrix element of field measurements. Besides, the distance which separates the magnetic sensors to the conductors is about a few millimeters in all these approaches. Two cases are studied, the first consists to reconstruct a small current I c in the presence of the earth's magnetic field (50 T), and in the second case, we add a disruptive conductor carrying a hard current I d (I d = 300 • I c ) at a distance of 3 cm to the system, the relative errors incurred for the reconstruction of the small current I c basing to the three preceding approaches are given in Table I.

The approach of field measurement leads to an inexpensive current sensor but it presents a poor magnetic background rejection. Whereas, by applying the field gradient measurement, the relative error obtained for the reconstruction of a small current is zero percent in the presence of the earth's magnetic field. This result is obvious because uniform magnetic disturbances are naturally subtracted. Then, we notice a better background rejection in the presence of the disruptive conductor I d (I d = 300 • I c ). One the side to these advantages, the gradient measurement solution is more expensive, and we will show in the following section, that the position of the directive magnetic sensor is very complex and it will affect on the computation of the relative error. However, the gradient of gradient measurements method improves significantly the rejection of an external magnetic disturbance, but the cost of the sensor therefore increases with using nine magnetic sensors and we note that the positioning of magnetic sensor is also critical. In the section three, we will detail the influence of misalignment of directive magnetic sensors on the computation of the relative error.

DIFFICULTIES OF IMPLEMENTING THE PRACTICAL APPLICATION FOR THE GRADIENT MEASUREMENT

When implementing the practical application for the approach of gradient measurements, many difficulties appear: the magnetic gradient should be measured along one single collinear axis; any offset introduces by the magnetic sensors will affect the computation of relative error and the rejection of disturbance fields is thus inadequate. By positioning a directive magnetic sensor at a distance d to the conductor crossing by the current to be reconstructed I c , to apply the approach of gradient measurement, we placed another directive magnetic sensor nearby the first magnetic sensor. The distance which separates the two magnetic sensors is . In practical application, an offset can be executed. If we suppose here, that the second magnetic sensor presents an offset in comparison to the initial setting (Fig. 3). The relative error due to this offset is not negligible as illustrate (Fig. 4). The relative error incurred for an offset of 1 is 2.3% in a cleaning magnetically environment and 2.74% in the presence of the disruptive conductor (I d = 300 • I c ). It's higher than the relative error required 2%. However, the subtracting of the two field measurements in 2D is given by (10):

B x = 0 I c 2 • - d + d ( 10
)
If c is the measured error committed on the subtracting magnetic field, the relative error for the reconstructing of I c is:

e = c • d • 1 + d (11)
As a conclusion, a thin distance between the magnetic field sensors will increase the percentage of the relative . Fig. 3. A misalignment of one directive magnetic sensor is executed (two magnetic field sensors are represented by two arrows, in dotted we present the initial position of the directive magnetic sensor; d is the distance which separate the first magnetic sensor to the conductor to be measured I c , and is the distance between the two magnetic sensors). error taking into count the measured error. Then, if we increase the distance which separates the magnetic field sensors, the approach of gradient measurement is not thus reach its aim; because to assure a good rejection of magnetic disturbance, the two magnetic field sensors should be placed closed to each other.

The Figure 5 shown the relative error incurred when the measured error c is equal to 1% and in a cleaning magnetically environment; more the distance between the magnetic sensors is higher the relative error decrease but it is lead to a poor magnetic background rejection. In the following sections, we will search a solution to realize the practical application for the approach of gradient measurement, able to subdue the disadvantages of the above sensor/conductors configuration and having a satisfactory cost.

IMPLEMENTATION OF GRADIENT MEASUREMENT WITH "ROUND-TRIP" CONDUCTORS

The approach of gradient measurements can be also applied by using other configurations. If we suppose, that the conductor is made in a round-trip configuration, and by positioning one field sensor above each conductor section (Fig. 6(a)), we can reconstruct the current flowing through the conductor by subtracting the two measured fields. The theoretical study for this configuration predicts the same results for rejection of magnetic disturbances. To get a relative error of 2% for reconstruction I c , the disruptive conductor I d (I d = 300 • I c ) should be placed at a distance higher than 13 cm. Figure 7 illustrates the model using six field detectors for three "round-trip" conductors to apply the gradient measurement method. The distance which separates each conductor sections is also a few millimeters. Once the geometry of the system is fixed, we can express the matrix for the approach of magnetic gradient measurement by subtracting two magnetic inductions for two field detectors collinearly along the same axis above the conductor. The new matrix form is given in (12).

I 1 I 2 I 3 S 1 S 2 3 5 S 5 S 6 S 3 S 4
Fig. 7. Three-phase sensors with configuration "round-trip".

    H 1 -H 2 H 3 -H 4 H 5 -H 6     =     M 11 M 12 M 13 M 21 M 22 M 23 M 31 M 32 M 33     •     I 1 I 2 I 3     (12)
Where M ij is the difference of two elements k ij in the matrix of field measurements. These coefficients k ij depend on the geometry of the conductors, the distance between the conductor and the sensors and the distance between the conductor sections.

In the presence of the earth's magnetic field, the relative error incurred for the reconstruction of a small currant is zero percent. This result assures the interest of the gradient measurement to reject the uniform magnetic disturbance. Otherwise, if we placed at 3 cm a disruptive conductor carrying a hard current value I d (I d = 300 • I c ), the relative error for I c decreases to 27% (From 1016% field measurement). The results show an improvement of the rejection of non-uniform magnetic disturbance.

This configuration subdues the problem of misalignment of magnetic sensor, because their positioning is very easy mechanically because each couple of magnetic field sensor will be mounted on the same layout. Then, we suppose that c is the measured error committed on the subtracting magnetic field given by the relation (13) where p is the distance between the conductor sections and d the distance that separate the magnetic sensors from the conductor sections.

Bx = 0 I c 2 • -2p 2 d p 2 + d 2 (13)
The relative error for the reconstruction of I c is more lower than the previous case, where the two magnetic sensors were placed closed to each other, and it is given by ( 14):

e = c • d • 1 2 + d 2 2p 2 (14)
As we conclude, a thin distance between the conductor sections will increase the percentage of the relative error taking into count the measured error, but with this configuration, the relative error is multiply by a factor 1/2 instead a factor 1 as in the preceding case. Then, this relative error will greatly decrease. Besides, if we increase the distance which separates the conductor sections, we will then increase the dimension of the current sensors what is conflicting with size reducing aim of this new technology. The Figure 8 shown the relative error incurred when the measured error c is equal to 1% and in a cleaning magnetically environment; more the distance between the conductor sections is higher the relative error decrease but the dimension of the current sensor rise against the preceding case.

We can conclude that the application of the approach of gradient measurement with "round-trip" configuration presents a good magnetic background rejection and subdue the disadvantages of the above sensor/conductors configuration having two magnetic sensors placed closed to each other nearby the conductor crossing by the current to be measured. Then, with round-trip conductor we can use a directive magnetic sensor without a problem of misalignment because it is placed on the same layout. Finally, with a distance of a few millimeters between the conductor sections, the relative error incurred is thus satisfactory. However, the addition of a simple magnetic shielding to the configuration "round-trip" can lead to the design of highly accurate and low cost current sensors, and in this case the material of shielding isn't saturate.

CONCLUSION

A method for non-contact current measurement based on magnetic gradient measurement is described. Several field probes in close proximity to a current carrying conductor are combined to measure the spatial derivative of the magnetic field, and which leads to a more accurate evaluation of the current. (We have demonstrated the advantages of this approach by the theoretical study for a one conductor case.) Furthermore, we have extended the methodology for the case of three phases; the numerical results have shown a much improved rejection of the homogeneous perturbations with the use of gradient measurements. In addition, in the presence of a disruptive conductor placed at 3 cm and carrying a hard current value (I d = 300 • I c ), we have reconstructed a small current I c with a lower relative error by applying the gradient approach. However, when implementing the practical application for the gradient measurement approach, many difficulties will be appearing as the misalignment of directive magnetic sensors and the influence of the distance between their on the computation of the relative error for reconstructing the currents. Next, we have proposed a solution able to subdue the upper difficulties; it is based on "round-trip" conductors and consists to laying the two magnetic sensors collinearly along the same axis above the conductor sections to apply gradient measurements. This configuration was showing a good magnetic background rejection and a satisfactory cost. Besides, with round-trip conductor we can use a directive magnetic sensor without a problem of misalignment because it is placed on the same layout. Finally, the addition of a simple magnetic shielding to the configuration "round-trip" can lead to the design of highly accurate and low cost current sensors.

Fig. 1 .

 1 Fig. 1. Exclusion areas for the disturbing conductor with < 0 02 for the three approaches. The field sensors, represented by an arrow, are positioned near to the conductor with the current to be measured I c ; I d is the hard current value in the disruptive conductor.

I 3 H 2 I 1 Fig. 2 .

 3212 Fig.2. n magnetic field sensors surrounding three conductors.

Fig. 4 .Fig. 5 .

 45 Fig.4. The relative error due to the misalignment of one magnetic field sensor.

Fig. 6 .

 6 Fig. 6. (a) Two magnetic field sensors represented by two arrows are positioned above the sections of "round-trip" conductor flowing by the current to be measured I c ; I d is the disruptive conductor carrying a hard current value. (b) Exclusion areas for disturbance conductor with < 0 02 for measurements.

Fig. 8 .

 8 Fig. 8. The relative errors function on the distance which separate the conductor sections.

Table I .

 I The relative error for the reconstructing of a small current I c in the presence of a homogeneous and non homogeneous magnetic field.

		Relative error (%)
		In the presence of	In the presence of a
		earth magnetic	disruptive conductor
		field (50 T)	(I d = 300 • I c )
	Inversion field	245	1016
	measurements
	Gradient measurements	0	55 6
	Gradient of gradient	0	4 05
	measurements

Acknowledgments:

The authors would like to thank "La METRO de Grenoble" who supported this work.


Chuang Wang

Benoît Tauzin
email: benoit.tauzin@univ-lyon1.fr

Thanh-Son Pham

Hrvoje Tkalčić

On the efficiency of P-wave coda autocorrelation in 1 recovering crustal structure: examples from dense 2 arrays in the eastern US 3

We investigate the efficiency of P-wave coda autocorrelation in recovering crustal 9 architecture 10 • We compare imaging from P-wave coda vertical autocorrelations and radial re-11 ceiver functions using dense profiles from US arrays 12 • Vertical autocorrelations clearly outperform radial receiver functions in presence 13 of a sedimentary basin

Introduction

Teleseismic earthquakes occur at distances of more than thousands of kilometers away from the recorders (typically epicentral distances larger than 30 • ). At these epicentral distances, the incoming teleseismic wave is approximately a plane wave with a relatively steep incidence angle. [START_REF] Burdick | Time functions appropriate for deep earthquakes[END_REF] noticed that vertically polarized teleseismic waveforms are dominated by the source, the near-source structure and the instrument effects. The near-receiver structure, due to a small incidence angle, influences mainly the compressional wave (P wave) response in the vertical seismograms.

Similarly, the shear wave (S wave) response of the structure beneath a receiver can be retrieved from the horizontal seismograms. [START_REF] Langston | Structure under mont Rainier, Washington, inferred from teleseismic body waves[END_REF] proposed a technique to determine the near-receiver structural response by deconvolving the vertical component of teleseismic seismograms from the horizontal ones. Through deconvolution, the complicated source and wavefield propagation effects away from the receiver are removed from the horizontal components, and the results are referred to as the radial receiver functions (RFs). Vinnik (1977) introduced slightly earlier the correlation-based approach to compute receiver functions.

The introduction of auto-and cross-correlations to seismology predated the development of the RF method. The autocorrelation of the transmission response of an acoustic plane-wave was used by [START_REF] Claerbout | Synthesis of a layered medium from its acoustic transmission response[END_REF] to extract the reflection response of the stratification beneath a receiver. Autocorrelation is a mathematical operation that measures the similarity of a waveform with itself in a single time series. A similar operation, crosscorrelation, involves two seismograms, and is used to build the radial RF with deconvolution [START_REF] Clayton | Source shape estimation and deconvolution of teleseismic bodywaves[END_REF]Vinnik, 1977;[START_REF] Ammon | The isolation of receiver effects from teleseismic P waveforms[END_REF]. In the following, we will interchangeably use the terminology "autocorrelation" and "vertical RF", because in practice the vertical RF can be built from the autocorrelation of the vertical component seismogram (Tauzin et al., 2019).

A few early studies used the principle of correlation for characterizing crustal structure (e.g., L. J. [START_REF] Burdick | Modeling crustal structure through the use of converted phases in teleseismic body-wave forms[END_REF] X.-Q. [START_REF] Li | Deconvolution of teleseismic body waves for enhancing structure beneath a seismometer array[END_REF][START_REF] Langston | The vertical component p-wave receiver function[END_REF]. In particular, C. [START_REF] Langston | The vertical component p-wave receiver function[END_REF] suggested that all three components of the RFs should be involved in inferring Earth-structure models, and demon-strated the reliability of the technique by means of using vertical-component P-wave (autocorrelation) RFs. More recent developments of the autocorrelation method make use of the vertically propagating component of the seismic ambient noise (e.g., Tibuleac & von Seggern, 2012;[START_REF] Gorbatov | Crustal properties from seismic station autocorrelograms[END_REF][START_REF] Kennett | Lithosphere-asthenosphere p-wave reflectivity across australia[END_REF]Kennett et al., 2015;Taylor et al., 2016;[START_REF] Saygin | Retrieval of the p wave reflectivity response from autocorrelation of seismic noise: Jakarta basin, indonesia[END_REF][START_REF] Becker | Crustal thickness across the transeuropean suture zone from ambient noise autocorrelations[END_REF]. The difference between the noise-based and event-based approaches is still under investigation [START_REF] Kennett | Lithosphere-asthenosphere p-wave reflectivity across australia[END_REF]Sun & Kennett, 2016;[START_REF] Romero | Mapping the basement of the ebro basin in spain with seismic ambient noise autocorrelations[END_REF][START_REF] Buffoni | Crustal structure beneath tierra del fuego, argentina, inferred from seismic p-wave receiver functions and ambient noise autocorrelations[END_REF]Tauzin et al., 2019). The choice of filtering have been shown to play an important role in the recovery and type of recovered signal from noise autocorrelations [START_REF] Helffrich | Comment on "Crustal thickness across the Trans-European Suture Zone from ambient noise autocorrelations" by G. Becker and B. Knapmeyer-Endrun[END_REF] with two main approaches. The first approach is inherited from standard radial RFs analysis and makes use of frequencies smaller than 1 Hz (Tibuleac & von Seggern, 2012;Taylor et al., 2016). This approach allows for the identification of individual reflected phases on autocorrelograms for subsequent picking and interpretation in terms of stratigraphy. The second approach makes use of the high-frequency component of the noise that involves an analysis within a narrow frequency band [START_REF] Gorbatov | Crustal properties from seismic station autocorrelograms[END_REF][START_REF] Kennett | Lithosphere-asthenosphere p-wave reflectivity across australia[END_REF][START_REF] Becker | Crustal thickness across the transeuropean suture zone from ambient noise autocorrelations[END_REF][START_REF] Romero | Mapping the basement of the ebro basin in spain with seismic ambient noise autocorrelations[END_REF]. In this second case, deterministic seismic arrivals are not anymore detected, and the interpretation is driven by changes in the reflectivity patterns in autocorrelations.

In sedimentary basins (SBs), radial RFs are contaminated by a strong resonance effect. As shown by e.g. Zheng et al. (2005), RFs from recording sites that are filled by sediments are dominated by reverberations in the sedimentary layers within the first few seconds after the direct P-arrival. Indeed, due to a strong, basal, impedance contrast, the low-velocity sedimentary layer (LVZ) generates quite strong P-wave and S-wave reverberations and conversions. The superimposed reflections and conversion at the base the sediments often dominate over the converted phases from deeper crustal discontinuities and mask corresponding information about Earth structure [START_REF] Frederiksen | Deriving crustal properties from the p coda without deconvolution: the southwestern superior province, north america[END_REF]. Thus, conventional methods based on the analyses of travel-times of converted waves can lead to erroneous interpretations.

Several methods have been developed for mapping the crustal structure in SBs. Those methods, designed to reduce SB reverberations, include downward continuation and wavefield decomposition [START_REF] Langston | Wave-field continuation and decomposition for passive seismic imaging under deep unconsolidated sediments[END_REF], autocorrelation with a resonance removal filter in the frequency domain (Yu et al., 2015), and radial RF method to fit the waveforms and invert for elastic structure of a SB [START_REF] Piana Agostinetti | Sedimentary basin investigation using receiver function: an east african rift case study[END_REF]. Wavefield downward continuation method [START_REF] Langston | Wave-field continuation and decomposition for passive seismic imaging under deep unconsolidated sediments[END_REF] requires a priori knowledge of Earth structure and modeling to extract wave conversions. This method has been applied to estimate sedimentary and crustal structures (e.g., Tao et al., 2014). Another recent study investigated sedimentary structure of basins in China from frequency dependent P-wave particle motion [START_REF] Li | Joint inversion of rayleigh wave phase velocity, particle motion, and teleseismic body wave data for sedimentary structures[END_REF]. [START_REF] Li | Joint inversion of rayleigh wave phase velocity, particle motion, and teleseismic body wave data for sedimentary structures[END_REF] reported a joint inversion of Rayleigh wave phase velocity, particle motion, and teleseismic body wave data for sedimentary structures.

Here, we investigate the performance of a new method for imaging the crustal architecture based on the spectral whitening and autocorrelation of the vertical component of teleseismic records [START_REF] Pham | On the feasibility and use of teleseismic p wave coda autocorrelation for mapping shallow seismic discontinuities[END_REF]Tauzin et al., 2019). We investigate whether this method is appropriate for characterizing the crustal structure in presence of a SB. Spectral whitening is the division in the frequency domain of the original spectrum by the running average of its amplitude. We show that this operation efficiently filters out undesirable resonance associated with SBs. In comparison to the inversion of radial RF waveforms, autocorrelation based on filtering is computationally efficient for achieving direct imaging of the basin and crustal structure.

In order to demonstrate our method with dense seismic arrays, we analyze broadband data from temporary seismic arrays with quasi-linear geometry in the eastern US.

Recently, two studies have used the principle of correlation to analyze the crustal structure in the US. [START_REF] Delph | Constraining crustal properties using receiver functions and the autocorrelation of earthquake-generated body waves[END_REF] constrained the Moho depth and crustal Vp/Vs ra-tio underneath 88 stations from the Transportable Array using vertical and radial component RFs. van IJsseldijk et al. ( 2019) used a dense seismic array in Minnesota and Wisconsin and a Global Phase Seismic Interferometry approach, which is based on autocorrelations of teleseismic P-waveforms, but recorded at small incidence angles.

The seismic arrays that we use for imaging the crust in eastern US run from the Quaternary deposits in the Mississippi embayment, through a region of thickened crust at the western border of the Appalachian, to Maine in the northeast. Our data are divided into three sub-arrays with various station density and operation time that sample different geological settings. We investigate in particular the efficiency of the method for imaging a laterally variable crust overlaid in places by SBs. We present the general features of the processing. Given the amount of data and the processing for source equalization, a data selection procedure is implemented, which is dedicated to improving the signal-to-noise ratio of specific reflected arrivals. We observed coherent signals along the three sub-arrays. A comparison with a reference crustal model suggests that we are able to unambiguously image the Moho and the internal stratification of the crust, in contrary to radial RFs that poorly perform in presence of a SB. We also find in the eastern US the Moho to be a transitional layer (up to 5 km thick) instead of a sharp boundary.

Methods and Data

An example of vertical component RF for a simple one-layer crust is shown in Fig- ure 1b, with labels for major converted and reflected phases P s, P P p, P P s and P Ss.

The paths are shown in Figure 1a. We also show in Figure 1b for comparison the theoretical response recorded on the radial component. In our particular implementation of the vertical receiver function method, the extracted signal corresponds to the sequence of reflections (Figure 1b) in the stratification underneath the seismometer (Figure 1a).

The main reflection being recorded on this vertical RF is the reflection from the top of the Moho discontinuity, PPmp, where "m" marks the Moho (Figure 1b). Its amplitude results from the interaction of the teleseismic wavefield with the free-surface, with a reflection coefficient close to -1, and with the Moho interface, with a positive reflection coefficient. The combined reflection response gives a negative amplitude for the Moho reflection on the vertical RF record (Figure 1b). This response is not exactly the reflection response as obtained within the frame of seismic interferometry [START_REF] Claerbout | Synthesis of a layered medium from its acoustic transmission response[END_REF]Wapenaar et al., 2010;[START_REF] Ruigrok | Global-phase seismic interferometry unveils pwave reflectivity below the himalayas and tibet[END_REF]Tauzin et al., 2019). Indeed, if we take the reflection response due to a virtual source located at the surface at the reflection point of the teleseismic P-wave, then the negative reflection at the free surface is removed through the principle of interferometry, and the remaining reflection at the Moho is positive. Our estimated reflection response has therefore an opposite (negative) sign with respect to the P-wave reflection response in seismic interferometry imaging applications.

We compute the vertical RFs using spectral whitening, followed by autocorrelation [START_REF] Pham | On the feasibility and use of teleseismic p wave coda autocorrelation for mapping shallow seismic discontinuities[END_REF]Tauzin et al., 2019). If Z(f ) is the spectrum of the vertical component, the vertical receiver function E z (f ) is obtained from

E z (f ) ≈ Z(f ) Z * (f ) Z (f ) , (1)
where the superscript * marks the complex conjugate and Z(f ) Z * (f ) is the autocorrelation written in the frequency domain. Z (f ) is the regularized vertical component power spectrum taken as an approximation for the teleseismic source function. Each frequency sample is obtained through smoothing over a window of width W so the regularized vertical spectrum is

Z (f ) = 1 W f +W/2 f -W/2 ||Z(ν)|| 2 d ν.
(2)

-4-©2020 American Geophysical Union. All rights reserved.

manuscript submitted to JGR: Solid Earth

In the post-processing stage, we apply a 4th-order Butterworth filter with a bandpass chosen depending on the target structure. The frequency content necessary to image the shallow structure is higher (0.4 -1 Hz band-pass) than for imaging the Moho (0.2 -1 Hz band-pass). We also use different values of the smoothing width W depending on the target interface (see [START_REF] Pham | Antarctic ice properties revealed from teleseismic p-wave coda autocorrelation[END_REF]. We discuss the choice of the W parameter in section 5.3. We use a constant value of W = 0.1 Hz for targeting the Moho. We set W to a larger value, 0.4 Hz, when targeting the shallow stratification.

To improve the signal-to-noise ratio, all the vertical component RFs for a single station are stacked to build a high-quality reflection response. As the RFs are obtained over a large range of epicentral distances, the travel-time of P-wave reflected phases follows a move-out due to the variable path-lengths of waves incoming on the stratification. For a single layer of thickness H and P-wave velocity v p overlying a homogeneous halfspace, the travel-time of the P P p mode with respect to direct P (Figure 1b) is expressed by

T P P p (H, v p , p) = 2H 1 v 2 p -p 2 , ( 3
)
where p is the slowness of the impinging P-wave. We use the generalization of this expression to horizontally stratified media to compute a time-to-depth relationship and convert the time-domain auto-correlations into the depth-domain (Tauzin et al., 2019). This relationship is mode-dependent, and leads to destructive interference of signals that do not follow the move-out of P-wave reflections P P p. Single station auto-correlations are then linearly stacked in the depth domain. Due to the fact that the PPms phase has a similar travel-time relation as PPmp, the PPms phase appears as a "second Moho" at deeper depth with weaker amplitude (see Z component in Figure 1b). We tested phaseweighted stacking [START_REF] Schimmel | Noise reduction and detection of weak, coherent signals through phase-weighted stacks[END_REF] but did not find a major improvement in the quality of crustal images. We apply bootstrap resampling [START_REF] Efron | Statistical Data Analysis in the Computer Age[END_REF] for estimating the variability of stacked waveforms.

Despite the fact that equation (1) can be modified to extract radial receiver functions (Tauzin et al., 2019), we computed classical radial receiver functions through waterlevel deconvolution [START_REF] Clayton | Source shape estimation and deconvolution of teleseismic bodywaves[END_REF]. The difference resides in the (quasi-)absence of the P P p mode in water-level receiver functions (Tauzin et al., 2019). In that case, the radial RF is given by

E r (f ) ≈ R(f ) Z * (f ) Z (f ) , (4)
where R(f ) is the radial component spectrum. The regularized vertical spectrum Z (f ) is the initial spectrum that is leveled to a minimum amplitude defined by the water-level c [START_REF] Clayton | Source shape estimation and deconvolution of teleseismic bodywaves[END_REF],

Z (f ) = max(||Z(f )|| 2 , c ||Z|| 2 max ). ( 5
)
We apply in this case a Gaussian filter [START_REF] Clayton | Source shape estimation and deconvolution of teleseismic bodywaves[END_REF][START_REF] Ammon | The isolation of receiver effects from teleseismic P waveforms[END_REF]. The radial RFs are then depth-converted using a time-to-depth relationship for the direct Ps conversion mode

T P s (H, v p , v s , p) = H 1 v 2 s -p 2 - 1 v 2 p -p 2 . ( 6
)
Radial RFs at a single station in the depth domain are then linearly stacked to build the conversion response.

These methods are applied to broadband stations located in the eastern US (Fig-

ure 1). This choice is motivated by the presence of seismic arrays with various station density and operation time that sample different geological settings. We investigate in particular the efficiency of the method for imaging a laterally variable crust underlying in places thick sedimentary basins (SB) such as in the Mississippi Embayment. Our data

-5-©2020 American Geophysical Union. All rights reserved. Vertical (black) and radial (red) responses of the ground below the receiver for a simple, onelayer crust. In the vetical RF, the PPmp (PPp) and PPms (PPs) phases both exists but PPmp is much stronger (Tauzin et al., 2019). c) Topography, together with seismic sub-arrays used in this study and the seismicity (black circles) from the USGS catalog. d) Main geological deposits color-coded as a function of age (P. G. [START_REF] Schruben | Geology of the conterminous united states at 1: 2,500,000 scale a digital representation of the 1974 pb king and hm beikman map[END_REF]. e) Seismic sub-arrays and crustal thickness taken from the globally compiled model CRUST1.0 [START_REF] Laske | Update on crust1. 0-a 1degree global model of earth's crust[END_REF]. f) Sedimentary cover taken from CRUST1.0. The arrays belong to the following networks: Sub-Array I (red triangles) TA, Transportable Array; Sub-Array II (blue triangles) X4, Gulf of Mexico Basin Opening; Sub-Array III (white hexagrams) TA, Transportable Array.

are separated in three sub-arrays (Figure 1). We use the global crustal model Crust1.0 [START_REF] Laske | Update on crust1. 0-a 1degree global model of earth's crust[END_REF] and the geology of the conterminous United States [START_REF] King | Geologic map of the united states[END_REF]; P. G. A. Schruben et al., 1994) to describe the a priori structure beneath these arrays (Figure 1). We benefit from the Transportable Array for constructing transects of the crustal structure over several thousand kilometers with an average inter-station The set of data used in this study is the same as used for classical radial RF studies. Metadata and data are provided by the Incorporated Research Institutions for Seismology (IRIS) data center (DMC) [START_REF] Ahern | Seed standard for the exchange of earthquake data reference manual format version 2.4[END_REF] and the Global Centroid Moment Tensor (GCMT) catalog [START_REF] Ekström | The global cmt project 2004-2010: Centroid-moment tensors for 13,017 earthquakes[END_REF]. We used three-component stations with the broadband channel BH, and selected events from the GCMT catalog with moment magnitudes larger than 5.5 and epicentral distances between 30 and 90 • to avoid the effect of triplications at mantle transition zone discontinuities and the core shadow zone.

The processing is implemented using Matlab R scientific programing language, with the routine available at http://perso.ens-lyon.fr/benoit.tauzin/data.htm.

The seismograms were extracted and preprocessed to remove asynchronous records on the three components, resampled, and projected to the desired vertical-radial-transverse component orientations. We selected data with seismogram signal-to-noise ratio larger than 1.7, the signal and noise being measured from a root mean square over windows following ([0 3.25] s) and preceding ([-2.5 -0.5] s) the theoretical P-wave arrival (at 0 s) in the iasp91 model [START_REF] Kennett | Travel times for global earthquake location and phase identification[END_REF]. The seismograms are pre-processed to start 20 s before the P-wave arrival with a duration of 80 s. The seismograms are then processed using equations ( 1) and (4), i.e. applying spectral whitening and autocorrelation to build vertical RFs and water-level deconvolution for radial RFs.

The time-to-depth relationships, i.e. equations ( 3) and ( 6), are computed based on the velocity structure in the Crust1.0 model [START_REF] Laske | Update on crust1. 0-a 1degree global model of earth's crust[END_REF] at corresponding station locations. Crust1.0 is not taken here as providing the ground truth, but serves for a comparison with the observed data as well as guiding the interpretation of the recovered reflectivity. Crust1.0 is a global model that provides crustal types, thickness, and elastic properties of the crystalline crust according to basement age and tectonic setting [START_REF] Laske | Update on crust1. 0-a 1degree global model of earth's crust[END_REF]. The Moho depth in particular is constrained by data from activesource seismic studies as well as radial RF analysis. So this model provides a good reference for calibrating our seismic images and discussing the deviation of the first order structure from this model. The model in each cell is described as an eight-layer profile, including three sedimentary layers and three crustal layers. The parameters including P-and S-wave velocity, v p , v s , density, ρ and thickness, H, are given for each layer. Crust1.0 is discretized over a 1×1 • pixel spatial grid [START_REF] Laske | Update on crust1. 0-a 1degree global model of earth's crust[END_REF]. We do not expect to recover exactly the same structure beneath stations that do not have contributed to the model. This is the case for array II, which is densely sampling the equivalent of ∼2.5 cells of the Crust1.0 model (Figure 1). In all our applications, we compare the result of processing the observed data to synthetic data computed based on this global crustal model.

We simulate the waveforms with the plane-wave matrix algorithm, "respknt" [START_REF] Kennett | Seismic wave propagation in stratified media[END_REF][START_REF] Randall | Efficient calculation of differential seismograms for lithospheric receiver functions[END_REF]. Synthetics waveforms are computed for the same acquisition ge-ometry as in the observed data and processed with the same strategy: filtering, autocorrelation, time-to-depth conversion, and stacking. This approach is applied in all our synthetic experiments (Figures 2,6, 7, 9e, 10, 11, S4-S7), except Figure 9a-9d where the slowness is 0.06 s/km.

For observed data, we found in a preliminary attempt that the selection based on a signal-to-noise criterion applied on seismograms alone is not sufficient to ensure good quality auto-correlation stacks. At frequencies investigated here, the noise on single autocorrelations is a combination of incoherent (ambient) noise and signal-generated noise.

The signal-generated noise is a result of the teleseismic source function: it includes the characteristics of the seismic source (slip history, radiation pattern, source spectrum) as well as the effects of wave-propagation near the source. In addition, interfering phases, such as PP or PcP, add complexity unrelated to the receiver-side structure. Stacking a large number of data is required to decrease the noise level to acceptable level (Tauzin et al., 2019). But because we use temporary arrays (Figure 1), the number of available data is relatively small, which limits the performance of the stacking operation. Our experiments showed that at stations with a small number of data, even an application of a stacking with weighting based on the instantaneous phase [START_REF] Schimmel | Noise reduction and detection of weak, coherent signals through phase-weighted stacks[END_REF] is not able to provide satisfying results. We therefore discarded all stations with less than 10 (for sub-arrays I and III) or 15 (sub-array II) available autocorrelations.

We also relied on a more discriminative approach for data selection based on the use of a "phasing depth diagram". We use trial depths of seismic interfaces (phasing depths) and select traces that constructively contribute to the stacked signal at this depth. We provide in Figure S1 an illustrative diagram of the procedure, and in Figures S2 andS3 two examples based on simple synthetic signals and field data, respectively. For a negative signal corresponding to an interface at X km depth, we first make a selection of the traces with a negative signal within ±5 km from X km depth. We then stack the selected traces with bootstrap resampling. Because the traces are depth-converted, a successful selection of traces contributing to build a signal at X km depth should contribute to align (and build after stacking) signals for discontinuities at other depths. We retain only the traces that contribute to these multiple depths. We show in Supplementary Material how we use such a diagram juxtaposing stacked depth domain traces as a function of trial depths (the phasing depth diagram) to select a distribution of good quality data for a given station (supplementary Figures S1 andS3). This selective stacking method is more efficient applied to vertical RFs than to regular radial RFs. The reason is that the method focuses on removing the effect of random noise, but not coherent noise components. In regular radial RFs, a coherent noise component dominates, which is the structure-generated signal due to reverberations in the shallow sedimentary cover. Supplementary tables S1-S3 provide for the three sub-arrays the final number of RFs obtained at each station.

Synthetic experiments

We present in Figure 2 the Crust1.0 structure and corresponding synthetic radial and vertical RFs for selected stations ordered by increasing thickness of the sedimentary cover in Crust1.0. Radial RFs have been computed through water-level deconvolution following equations ( 4)-(5) (Figure 2b). Vertical RFs are calculated using autocorrelation and spectral whitening following equation ( 1)-(2) (Figure 2c). Generally, the vertical RFs provide good predictions for the Moho depth for all selected stations for all the velocity models. However, the radial RFs are more complicated. For a crustal cover without sedimentary basin (P56A), the 'Moho pulse' appears at the true depth and is reliable. For a thin sedimentary cover (QUA2), oscillations due to reverberations are very strong at shallow depth, while the 'Moho pulse' can still be distinguished in the radial RF. The converted pulses are followed by short duration and low-amplitude oscillations -8-©2020 American Geophysical Union. All rights reserved.

Results

Results for the three sub-arrays

We show in supplementary Figures S5,S6, and S7 examples of vertical RF stacks for three selected Transportable Array stations from sub-array I. For these stations, and despite a two-year duration of acquisition, there is no apparent loss in quality of the stacked autocorrelograms with respect to a permanent station such as HRV (Figure S4). The data selection slightly improves the lateral continuity of signals when building a profile from stations with variable numbers of earthquake records (Figure S8). For 71 stations in sub-array I, 10 stations have a set of less than 15 RFs, which are discarded, and among the remaining 61 stations, 43 stations have no more than 50 waveforms. The application of bootstrap resampling allows to estimate the variability of stacked waveforms and demonstrates that we observed robust signals (supplementary Figure S9).

Overall, the vertical RFs are relatively simple with 2, at maximum 3, major pulses (Figure 3). This is in contrast with the radial RFs, which show significantly more complexity (see section 5.2 in discussion; Figure 7). Chains of pulses are consistently observed along the three profiles at ∼30 km and ∼40 km depth (Figure 3 A priori geological and seismological information show that major holocene to cenozoic sedimentary basin structures are sampled by sub-array I and II (Figure 1d,f). Yet, we do not observe significant lateral variations of signal quality when the arrays cross such structures (Figure 3a). We will show in section 5.2 that this is a clear advantage of the vertical RFs over the traditional radial RF method. However, there are differences (Figure 4): (i) there exist minima and maxima in the vertical RFs that are not associated with discontinuities as described by Crust1.0, in particular, in the lithospheric mantle; (ii) many small amplitude arrivals might correspond to features present in the Crust1.0 model but are slightly offset, (iii) along the 250 km extent of the linear sub-array III, the data reveal variations of Moho topography (Figure 4b) that are not present in the model, (iv) at depths < 10 km, the response loses its interpretability as it appears dominated by signal-processing artifacts.

For further comparison, we computed the synthetic vertical RFs for the Crust1.0 model with the same acquisition geometry as observed in the data (Figure 6). This synthetic experiment demonstrates that vertical RFs are fairly simple and, in theory, most of the predicted pulses can be deterministically attributed to interfaces in the input model, although some secondary arrivals at crustal depth and in the uppermost mantle can be -12-©2020 American Geophysical Union. All rights reserved. attributed to "peg-leg reverberations" or cross-mode contamination (PPms). In addition, this synthetic experiment shows that the method, with this choice of parameters (f min = 0.2 Hz and W = 0.1 Hz), suffers from a loss in resolution at shallow depth (< 10 km). This is the result of the pronounced side-lobe of the peak at zero-lag in the autocorrelations.

We can draw three important conclusions from the comparison of Figures 4 and 6. Firstly, observed vertical autocorrelations are more complex than predicted based on Crust1.0, suggesting that there exist at crustal level small-scale heterogeneities and/or 3D effects of wave-propagation that are not captured using the reference model. Secondly, our method successfully recovers the crustal information without loss of signal-to-noise ratio or increased complexity in regions with thick sedimentary cover. This latter conclusion is supported by both the synthetic and observed data in the southwest portion of sub-array I and eastern portion of sub-array II (Figures 4 and6). Finally, in the observed data, P-wave reflections from the top of the Moho discontinuity are consistently observed with lower amplitudes than predicted by the Crust1.0 model (Figures 4 and6). In contrast, reflected pulses associated with intra-crustal discontinuities have a similar or slightly higher amplitude than predicted by Crust1.0. Supplementary Fig. S10 shows a negative correlation (-0.58) between our observed crustal thickness and the thickness of the sedimentary cover in the Crust1.0 model. This anti-correlation is predicted for a crust in isostatic equilibrium. The fit of a theoretical model to our best fitting line is obtained for a sediment density of 2.15 g/cm3, which is within the continental average range of densities for soft (2.07 g/cm3) and hard (2.38 g/cm3) sediments in the Crust1.0 model [START_REF] Laske | Update on crust1. 0-a 1degree global model of earth's crust[END_REF].

-13-©2020 American Geophysical Union. All rights reserved.

manuscript submitted to JGR: Solid Earth 5 Discussion

On the interpretation of autocorrelations

As presented in our results section, we unambiguously identify two chains of pulses near 30 and 40 km depth. From synthetic experiments and using of independent information (Crust1.0 model), we demonstrated that these chains of pulses are related to the reflections from the upper side of the Moho discontinuity, namely PPmp and PPms, and from the top of the lower crustal layer. Due to the limited number of data and intrinsic characteristics of the Moho boundary (see section 5.4), it can, however, be challenging to unambiguously pick the PPmp phase on autocorrelations. This is not a unique problem to our particular dataset, as most other prior studies based on autocorrelations suffer from similar interpretation problems (cross-mode contaminations and narrow frequency band limiting the resolution of velocity gradients) (e.g., Tibuleac & von Seggern, 2012;[START_REF] Gorbatov | Crustal properties from seismic station autocorrelograms[END_REF]Sun & Kennett, 2016;Taylor et al., 2016;[START_REF] Saygin | Retrieval of the p wave reflectivity response from autocorrelation of seismic noise: Jakarta basin, indonesia[END_REF][START_REF] Becker | Crustal thickness across the transeuropean suture zone from ambient noise autocorrelations[END_REF]. In addition, we show in next section that the vertical RFs behave better than radial RFs in the presence of a sedimentary cover. We therefore cast our analysis as a reference for future interpretations based on this technique.

Tauzin et al. ( 2019) showed that the lateral sampling of the PPp mode increases with the depth of investigation, and the difference with typical Ps conversion modes used in the radial RF imaging is of the order of 5 km at the depth of the Moho. This means that, in the single-station stacking, the cone of sensitivity of the P-wave reflections surrounding the stations will have a 10 to 20 km diameter (this distance is slowness-dependent) while it will have only a 5 to 10 km diameter for the Ps conversions. As such, stationbased gathering will be slightly more vulnerable to the presence of strong lateral variations in the Moho depth (crustal thickness) than in the radial RF imaging.

The effect will be a slight decrease in amplitude of the reflected phase due to outof-phase stacking. In addition to lateral variations of the Moho, crustal heterogeneity and surface topography both contribute to the PPmp-P differential travel time. The maximum topographic difference in our study is about 1 km. We estimate the delays introduced by a 1-km topography of the Earth surface and at the Moho discontinuity to be 0.8 s and 0.3 s respectively. Such delays are significant, but smaller than the width of the reflected pulse (1s). Consequently, if topographic variations over the area surrounding the station may contribute to the decrease in amplitude of the PPmp reflected pulse, they are unlikely to explain the whole difference with Crust1.0. This conclusion is also supported by the observation that the PPmp amplitudes are remarkably consistent across all of our images, and appear uncorrelated with local gradients in surface topography.

An additional strong argument for a negligible effect of surface topography is the observation of the reflections from the top of the lower crust with amplitudes similar as predicted by model Crust1.0. We conclude that the difference between the observed and the Crust1.0 Moho amplitudes must therefore have another origin (see section 5.4).

In supplementary Figure S11, we present an image of the crustal structure along the sub-array II based on the common reflection point stacking method (Tauzin et al., 2019), accounting for the geometry of illumination of the structure by the PPp mode.

Sub-array II is better positioned for performing such an experiment as the inter-station distance is the smallest among the three subarrays considered here. The image remains consistent with the observed structure from the single-station stacking, although the portions of the image with poor ray coverage suffer from a low signal-to-noise ratio.

Radial receiver functions

We present in Figure 7 7a-d) with unreliable converted pulses associated with the Moho discontinuity. This behavior is easily explained by the synthetic experiment shown in section 3. In conclusion, even in absence of noise of natural origin, the radial RFs are less reliable tools than the vertical RFs when low-velocity sedimentary layers exist near the surface. In comparison, the vertical RFs are simpler to interpret, although it is fair to say that they are not void of artifacts related to these shallow reverberations. We discuss the effect of reverberations and sedimentary basin structure in section 5.3, together with the impact of playing with the width of the smoothing window W for spectral whitening. We show in particular that a well-chosen value of W significantly improves the readability of the crust because spectral whitening can act as a resonance removal filter.

Filtering out basin resonance through spectral whitening

In sedimentary basins, a strong resonance effect is recorded on the vertical component as P-wave reverberates multiple times in the low-velocity sedimentary cover. Figure 8 shows the sequence of P-wave reverberations in a simple sedimentary structure, with a single layer. The sedimentary layer has lower velocities than the half-space beneath it. Teleseismic P waves reach the basement of the sedimentary layer and transmits as first-order P waves. The reflection of P waves at the free surface is a total reflection. Due to the significant impedance contrast at the base of the sedimentary layer, the free surface reflected P-waves reverberate multiple times inside the sedimentary layer.

These reverberations contribute to a strong resonance effect, which is recorded by the receiver on the Earth's surface mainly on the vertical component. In the time domain, these reverberations correspond to strong amplitude pulses that contaminate the seismograms. A remark is that on the vertical component, interference effects from converted S-waves are small because of the extremely low-incidence angle of P-waves due to the slow velocities in the SB layer. To understand the reverberations caused by SBs, we use an example with a sedimentary layer over a half space bedrock. The power spectrum of an impulse wavefront at normal incidence through a sedimentary layer over a half space bedrock can be expressed as,

|Z(f )| = 1 -r 1 + r • 1 (1 + (re -2πf Q ) 2 -2re -2πf Q cos 2πf τ ) 1 2 . ( 7
)
where details can be referred to Appendix A in [START_REF] Pham | Antarctic ice properties revealed from teleseismic p-wave coda autocorrelation[END_REF] for the full derivation of Equation 7 using the matrix propagation method [START_REF] Claerbout | Synthesis of a layered medium from its acoustic transmission response[END_REF][START_REF] Kennett | Seismic wave propagation in stratified media[END_REF].

In Equation 7, r = (ρ s v s -ρ c v c )/(ρ s v s + ρ c v c ) is the reflection coefficient at the sedimentary basin-crust interface in downward direction. ρ and v are density and wave speed of the media, where subscript s refers to the sedimentary basin and c to the crustal rock. τ = 2h/v is the two-way travel time in the layer. Q is the attenuation factor of in the sediment layer.

The spectral response in Equation 7 is a function in the frequency domain, which is characterized by peaks that are equally spaced at resonance frequencies,

f n = n + 1 2 • 1 τ . (8)
Despite the fact that Equation 7 is derived for normal incidence, it represents a good approximation for slight oblique incidence (see Figure S12 in the supplement). This shows that the contribution of converted S wave on the vertical component in slightly oblique incidence is negligible.

In equation ( 2), the frequency-dependent spectral-whitening operation allows balancing the frequency content of the RFs [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF][START_REF] Pham | On the feasibility and use of teleseismic p wave coda autocorrelation for mapping shallow seismic discontinuities[END_REF]Tauzin et al., 2019), and also removing specific spectrum-maxima related to the resonance of the layered structure beneath a receiver. To illustrate this, we proceed with a numerical experiment with simple one-dimensional models (Figure 9). The first model has a one-layer crust (no SB). The second model has a two-layer crust including a layer near the surface with low-velocities mimicking the effect of a sedimentary cover (SB model).

This layer is 2 km thick. The Moho is located at 30 km depth in both models. The elastic parameters chosen for the SB layer and the contrast in impedance across its base are consistent with the structure given for station QUA2 in Crust1.0.

We simulate vertical-component waveforms with the same plane-wave matrix algorithm [START_REF] Kennett | Seismic wave propagation in stratified media[END_REF][START_REF] Randall | Efficient calculation of differential seismograms for lithospheric receiver functions[END_REF]. Figure 9b shows the results of the numerical experiment in the time-domain for vertical component waveforms. In addition to the peak at zero-lag, the black trace (no SB model) shows the pulses of the P-wave reflection at the Moho (PPmp). For the SB model (red trace), the large contrast of elastic properties at the base of the SB introduces a strong resonance effect which contaminates the seismic response. This resonance corresponds to the ringing due to P-wave reverberations in the shallow sedimentary layer (Equation ( 7)). This makes it difficult to distinguish the reflection at the base of the crust (Figure 9b). In the frequency domain (Fig- ure 9c), the fundamental-mode resonance on the red spectrum appears as a broad and strong amplitude peak centered at ∼0.3 Hz, consistent with the 2-km thick, near-surface layer according to equation (3).

Spectral whitening offers a way to turn off such a resonance. The first step consists in smoothing the vertical component spectrum running an average over the pre-defined width W . We use W = 0.1 Hz in this example (Figure 9c). The smoothed spectrum (green) only preserves the broad peaks due to resonance in the shallow sedimentary layer.

The deconvolution in the frequency domain (Equation ( 1)) gives a spectrum close to the model without SB (Figure 9d), and the exponentially-decaying ringing is partially removed from the corresponding time-domain trace (Figure 9b). The whitened SB trace can now be interpreted, in particular to find the phase reflected at the Moho discontinuity. A residual ringing persists because the spectral whitening does not completely filter out the effect of the shallow SB on the vertical component spectrum. The same authors used the H-κ grid searching technique to study the Moho depth and average Vp/Vs ratio. They applied H-β method to remove the reverberations caused by unconsolidated sediments. Yang and Niu (2019) studied sedimentary structure of basins in China from frequency dependent P-wave particle motion. [START_REF] Li | Joint inversion of rayleigh wave phase velocity, particle motion, and teleseismic body wave data for sedimentary structures[END_REF] successfully isolated the Moho P-to-S conversion from sediment reverberations with a wavefield-downwardcontinuation technique. They reported a joint inversion that can retrieve the input sedimentary models. A comprehensive comparison between these methods and the resonance removal filter described here would be interesting but is out of the scope of this study.

The Moho as a transitional layer

In high-frequency (> 10 Hz) active seismic reflection experiments, it is common to observe a laterally variable crustal reflectivity between geological terrances, but also a vertical contrast in reflectivity at the transition between the crust and upper mantle (e.g., Kennett & Saygin, 2015). The Mohorovičić discontinuity is regularly identified as the base of the lower crust, but is often itself not marked by a strong reflectivity. The upper mantle shows no distinct reflections probably because the vertical scale of the structure is longer than the seismic wavelength used in active reflection surveys [START_REF] Kennett | Multiscale seismic heterogeneity in the continental lithosphere[END_REF]. At the same time, seismic surveys based on radial RFs with lower frequency content (< 1 Hz) point at the Moho being the most reflective interface.

The Moho is a first-order discontinuity in Crust1.0 with contrasts of ∼16% in v p and ∼13% in v s and density. Here, we find that amplitudes of low-frequency (0.2-1.0 Hz) P-waves reflections on top of the Moho, measured using autocorrelations (Figure 4), are weaker than predicted by the Crust1.0 model (Figure 6). We discussed in section 5.1 the effects of surface and boundary topography, pointing that these are unlikely explanations for such weak amplitudes. We here consider three alternative scenarios: (i) the elastic impedance contrast at the boundary is smaller than predicted by Crust1.0, (ii) there is an effect of anelasticity that is not accounted for in the seismic modeling, and/or (iii)

the Moho is a transitional layer rather than a sharp boundary. The first scenario is un- 2005; Thybo & Artemieva, 2013). We demonstrate here that a gradual transition could jointly explain our observations of amplitudes for topside P-wave reflections from vertical autocorrelations and P-to-S transmissions on radial RFs. P-wave reflections and P-to-S transmissions have a different sensitivity to velocity gradients (e.g. [START_REF] Bostock | Seismic waves converted from velocity gradient anomalies in the Earth's upper mantle[END_REF]Tauzin et al., 2019). We illustrate this by plotting the transmission and reflection coefficients T ps and Rpp as a function of the vertical extent h of a transitional layer adopting the contrasts in elastic properties for the Moho from Crust1.0 (Figure 10a). For this computation, we decompose the gradient in thin homogeneous layers and use a planewave matrix algorithm [START_REF] Chapman | Yet another elastic plane-wave, layer-matrix algorithm[END_REF][START_REF] Ma | A note on the equivalence of three major propagator algorithms for computational stability and efficiency[END_REF]Tauzin et al., 2019). We use a reference frequency of 0.6 Hz, which is the center of the frequency band of analysis for long-period vertical RFs. We also show the coefficient dependence on various slowness.

Generally, the transmission coefficients T ps decrease smoothly from a sharp interface to a ∼10 km-thick interface. In contrast, the reflection coefficients Rpp decrease abruptly and non-linearly with h from ∼10% to 0% and reach at 0.6 Hz a minimum for a ∼ 8 kmthick layer. This means that a transitional layer tends to affect the reflectivity of P-wave reflections more than the reflectivity from direct P-to-S conversions [START_REF] Bostock | Seismic waves converted from velocity gradient anomalies in the Earth's upper mantle[END_REF]Tauzin et al., 2019).

SW

Conclusion

We investigated the efficiency of an autocorrelation technique to extract vertical component receiver functions from teleseismic earthquake data and create images of the crust beneath dense seismic arrays. Using synthetic and field-data examples, we compared the performances of both vertical and radial RFs. We examined in particular the contribution of reverberations in shallow sedimentary (SB) layers to the images formed from both data types. Vertical RFs are simpler to interpret than radial RFs, although they are not void of complexities associated with shallow reverberations. We show that the operation of spectral whitening is critical for obtaining good quality images beneath SBs, as it acts as a resonance removal filter. We processed broadband data from three temporary USArray sub-arrays with quasi-linear geometry in the eastern US. By identifying pulses of reflected P-waves on the vertical autocorrelograms, we were able to identify the Moho and a laterally-variable crust structure beneath these subarrays. We observed coherent signals along the three sub-arrays despite the presence of SBs. Compared with the crustal model Crust1.0, our analysis suggests the existence of additional complexity related to fine-scale layering and/or 3D structure. We developed a method to constrain the thickness of the Moho boundary from an amplitude analysis of reflected and converted phases on vertical and radial component RFs. This analysis shows that the Moho beneath eastern US is a transitional layer extending over a depth interval of up to 5 km. The proposed method is highly complementary to other receiver-based methods and has the potential to advance our understanding of the crustal architecture with or without SBs.

  Figure 1. a) P-wave reverberations from teleseismic records. The seismic stations (black inverted triangle) sits on top of the crust (orange) above the mantle (brown). The thickness of the crust is denoted as H. P-wave is marked by solid lines and S-wave by dashed lines. The phases most relevant in this study are shown as Ps (blue), PPp (red), PPs (green), and PSs (black). b)

  spacing of 70 km (sub-arrays I and III). We distinguish Sub-Array I from III because they sample different geological structure. Sub-array I samples the thick sedimentary cover (up to 11 km) and thin crust (25-30 km) in the Mississippi embayment at southwest, and the western border of the Appalachian in the middle with a thicker crust (up to 50 km) and thinner sedimentary cover (3 km-thick) whereas sub-array III in Maine samples a standard crust (30 km thick) with no sedimentary cover. Finally, the Gulf of Mexico Basin Opening experiment (array II) is a linear array with 15 km inter-station spacing sampling sharp lateral changes of crustal and sedimentary thickness. Over 250 km, the model Crust1.0 indicates that the crust thins from 50 km at northwest to 35 km at southeast, and the sedimentary cover increases in thickness from 0 to 11 km. The 15-km inter-station spacing of sub-array II allows testing another imaging principle, Common Reflection Point (CRP) migration (Tauzin et al., 2019), which unfolds the information on the structure by projecting amplitude on the RFs onto the final leg of the path of reflected waves underneath the receivers (supplementary Figure S11).

Figure 2 .Figure 3 .

 23 Figure 2. Reference P-wave velocity profiles (a) and corresponding synthetic radial (b) and vertical receiver functions (c) for 10 selected stations with increasing thickness of the sedimentary cover according to Crust1.0. The red lines indicates the Moho depth in Crust1.0. The vertical RFs are processed with W = 0.1 Hz.

  ), likely associated with the deepest crustal layering. Yet, it can be relatively complicated to interpret vertical RFs in terms of a deterministic structure. Multiple arrivals with similar amplitudes are recovered at crustal levels, making the identification of the Mohorovičić discontinuity ambiguous. In some cases, arrivals are detected at lithospheric mantle depth, for example at the southwestern portion of sub-array III (Figure3c), where the Moho is not expected -10-©2020 American Geophysical Union. All rights reserved.manuscript submitted to JGR: Solid Earth to lie deeper than 35 km according to Crust1.0 (Figure1e). Although we cannot take the structure provided by Crust1.0 as the ground truth, we can interpret this deep signal as a fictitious interface that results from the arrival of multiple reverberations with cross-mode coupling, PPs (see Figure1a, b and synthetic experiment based on Crust1.0 in section 4.2).

Figure 4 .Figure 5 .

 45 Figure 4. Superimposition of observed P-wave auto-correlation signals (black traces) to the P-wave velocity structure (background colors) in the reference crustal model Crust1.0 (Laske et al., 2013). The white circles are associated with the local seismicity (USGS catalog; see Figure 1c). The layers are sed. (sedimentary layers), u.c. (upper crust), m.c. (middle crust), l.c.(lower crust), and m. (mantle).

  Figure 6. Superposition of synthetic P-wave auto-correlation signals (black traces) to the P-wave velocity structure (background colors) in the global crustal model Crust1.0 (Laske et al., 2013). The light purple and blue horizons are same as Figure 5.

Figure 7 .

 7 Figure 7. Observed (a, c, and e) and synthetic (b, d, and f) radial receiver functions (black traces) superimposed to the P-wave velocity structure (background colors) in the global crustal model Crust1.0 (Laske et al., 2013). Panels the others are synthetics. No post-processing filtering is applied.

Figure 8 .

 8 Figure 8. Sequence of reverberations in a sedimentary layer. The seismic station (black inverted triangle) sits on a sedimentary layer (orange) above the solid crust (brown). The direct P wave reaches the discontinuity between the sedimentary layer and the solid crust, and generates a first-order P wave. Due to the significant impedance contrast, the first-order P waves reverberate multiple times inside of the sedimentary layer. The thickness of the sedimentary layer is denoted as h. θP c and θP s are the incident angles of P waves in the crust and sediment layer respectively.

Figure 9 .

 9 Figure 9. Numerical experiment illustrating turning off the resonance associated with a shallow sedimentary layer. (a) Crustal velocity models with a "normal" crust (no SB) and a sedimentary basin (SB model). (b) Simulated vertical responses aligned by P-wave arrival. The black seismogram is the vertical response for the normal crust (no SB). The red one is for the sedimentary basin (SB model). The blue one is the result of spectral whitening for the SB model (SB white). Major identified phases are the P-wave reflections (PPmp) and the P-to-S reflection (PPms) above the Moho. The slowness p = 0.06 s/km. (c) Spectra Z(f ) of the vertical component responses. The spectra in black and red represent the responses ZZ * for the normal crust (no SB) and sedimentary basin (SB model) respectively. The spectrum in green represents the smoothed response Z (f ) for the SB model. (d) The spectrum after spectral whitening for the SB model (blue), compared with the spectrum for the "normal crust" model (black). (e) Synthetic vertical autocorrelations for the same 10 stations as in Figure 2 computed with W =0.4 Hz. Combined with Figure 2c, this figure shows filtering-on or -off the resonance in the shallow sedimentary cover with different value of smoothing parameter W .

Figures

  Figures 9e andFigure2c compare synthetic vertical autocorrelations for the var-

Figure 10 .

 10 Figure 10. Sensitivity of PPmp reflections and Pms transmissions to a gradual Moho layer. (a) Variation of the topside P-wave reflection coefficient Rpp (solid lines) and bottom-up transmission coefficient Tps (dashed lines) with the thickness of the velocity gradient at the Moho. Coefficients are calculated for a P plane wave with different slowness: p = 0.04 s/km (red), p = 0.065 s/km (black), and p = 0.09 s/km (blue). Here, the elastic properties across the Moho is defined by Crust1.0 for station KVTX in sub-array I. The value of coefficient is taken in the middle of the frequency band used in an analysis of vertical RFs (0.2 -1.0 Hz). (b) Histogram of measured amplitude of reflections (AP P p) for sub-arrays I and III. (c) Histogram of measured amplitude of conversions (AP s) for sub-arrays I and III. (d) Measured amplitudes of reflections (AP P p) as a function of amplitudes of transmissions (AP s) for observed data of sub-arrays I (blue circles) and III (red triangles). (e) AP P p v.s. AP s from synthetic (rectangle and diamond points) and observed data (grey circles and triangles) with Moho thickness h = 0 km and h = 6 km. Data belongs to sub-arrays I and III. (f) Theoretical curves of AP P p as a function of AP s, with points indicating the Moho thickness (h = 0, 5, 10, 15 km). The grey points are extracted from observed data that is same as (c).

Figure 11 .

 11 Figure 11. Superposition of synthetic P-wave auto-correlation signals (black traces) to the P-wave velocity structure (background colors) in the global crustal model Crust1.0 (Laske et al., 2013) modified using a transitional Moho layer 6 km thick.

  manuscript submitted to JGR: Solid Earth thick gradient, from which we reported the measurements shown in Figure10d. The predicted Moho amplitudes appear to be in better agreement with the observed data than the predictions from a first-order discontinuity (Figure6). We note that the Moho transition can also be reduced to a value thinner than 5 km if we consider a combination of possible mechanisms that can reduce the Moho reflection amplitudes, such as crustal transmission losses, 3D scattering, and effects of inelastic attenuation. However, we would have to explain why such mechanisms do not affect similarly affect the reflections from the top of the lower crust.

©2020 American Geophysical Union. All rights reserved.manuscript submitted to JGR: Solid Earth

Acknowledgments

The IRIS Data Management Center was used for access to waveform data (available from http://service.iris.edu/fdsnws/dataselect/1/). IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681.

We acknowledge Crust1.0 model (https://igppweb.ucsd.edu/ ~gabi/crust1.html).

We acknowledge the comments from reviewer Elmer Ruigrok, Associate Editor Sebastien Chevrot, and an anonymous reviewer. B.T. has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 793824. The Survey; For sale by USGS Map Distribution,. Schruben, P. G. A., Bawiec, R. E., King, W. J., Beikman, P. B., & Helen, M. (1994).

Geology of the conterminous united states at 1: 2,500,000 scale-a

-24-©2020 American Geophysical Union. All rights reserved.