Li-Yu Kan

Sébastien Chevrot
email: sebastien.chevrot@get.omp.eu

Vadim Monteiller

Dehydration of the Subducting Juan de Fuca Plate and Fluid Pathways Revealed by Full Waveform Inversion of Teleseismic P and SH Waves in Central Oregon

The dehydration of subducting slabs expels a massive amount of water into the forearc and backarc mantle which is responsible for the serpentinization of the mantle wedge, as well as the production of melt and arc magmatism. These processes are expected to have characteristic signatures in density and seismic velocities models, which remain largely elusive to date due to the limited spatial resolution of classical passive tomographic approaches. Here we present a tomographic model of density, V P , V S , and V P /V S beneath central Oregon, obtained by inverting complete teleseismic P and SH waveforms recorded by the CASC93 temporary experiment. The final model shows an east-dipping low-velocity layer less than 10 km thick that can be associated with the fluid-saturated Juan de Fuca oceanic crust. The distribution of tremors at the surface closely coincides with the horizontal extent of this low-velocity layer. Below 40 km depth, seismic velocities and density increase progressively to the values of a typical mantle. This transitional domain corresponds to the eclogitization of the oceanic crust. Silica-saturated fluids released by pore collapse migrate upward, producing serpentinization reactions in the forearc mantle that lower the density and seismic velocities. The very low V P / V S ratio documented in the Cascadia forearc crust is evidence that these silica-saturated fluids reach the crust, where they produce extensive quartz mineralization. At greater depth, a low seismic velocity and high V P /V S ratio anomaly provides evidence for partial melting at around 75 km depth beneath the volcanic arc.

Plain Language Summary

The subduction of the hydrated oceanic lithosphere injects large amounts of water into the Earth's interior. However, owing to the limited spatial resolution of classical tomographic approaches, the fluid pathways and related metamorphic reactions remain elusive. In this study, we image the Cascadia subduction beneath central Oregon by inverting complete P and SH teleseismic waveforms. The resulting tomographic models of density and seismic velocities reveal the eclogitization reactions in the subducting Juan de Fuca oceanic crust and the associated release of fluids into the forearc mantle. These silica-saturated fluids migrate upward, producing serpentinization of the forearc mantle. When these fluids reach the crust, the temperature decrease produces quartz mineralization, as evidenced by the very low V P /V S ratio observed in the Cascadia forearc crust. At greater depth, we find evidence of partial melting around 75 km depth below the volcanic arc, resulting from dehydration of the slab that lowers the solidus temperature in the lithospheric mantle of the upper plate. KAN ET AL.

Introduction

Subduction of hydrated oceanic lithosphere introduces massive amounts of water into the deep Earth (e.g., S. A. [START_REF] Peacock | Fluid processes in subduction zones[END_REF]. This water controls the mechanical properties of the subduction interface and the depth distribution of earthquakes and episodic tremor and slip (S. M. [START_REF] Peacock | Hydrous minerals in the mantle wedge and the maximum depth of subduction thrust earthquakes[END_REF]. The temperature-dependent transition from blueschist to eclogite facies governs the dehydration of the subducting plate (S. M. [START_REF] Peacock | Seismic consequences of warm versus cool subduction metamorphism: Examples from southwest and northeast Japan[END_REF]. For example, for the cold and old Pacific plate subducting beneath NE Japan, this transition is around 110 km depth whereas for the hot and young Philippine Sea plate subducting beneath SW Japan, the downgoing oceanic crust enters the eclogite facies at around 50 km depth (S. M. [START_REF] Peacock | Seismic consequences of warm versus cool subduction metamorphism: Examples from southwest and northeast Japan[END_REF]. These estimates are in good agreement with the seismologically constrained depth extent of a low seismic velocity layer often interpreted as the fluid-saturated oceanic crust (e.g., M. [START_REF] Bostock | The Moho in subduction zones[END_REF]. The eclogitization of the downgoing oceanic crust and the related collapse of porosity release large amounts of fluids in the upper plate. These silica-saturated fluids will percolate through and transforms the forearc mantle into weak assemblages of serpentine, talc, and brucite (S. M. [START_REF] Peacock | Hydrous minerals in the mantle wedge and the maximum depth of subduction thrust earthquakes[END_REF]. In warm subduction zones, the downdip limit of subduction earthquakes is thus controlled by the depth at which the subducting slab intersects the forearc mantle (S. M. [START_REF] Peacock | Hydrous minerals in the mantle wedge and the maximum depth of subduction thrust earthquakes[END_REF][START_REF] Oleskevich | The updip and downdip limits to great subduction earthquakes: Thermal and structural models of Cascadia, south Alaska, SW Japan, and Chile[END_REF][START_REF] Hyndman | Downdip landward limit of Cascadia great earthquake rupture[END_REF][START_REF] Hyndman | Cascadia subducting plate fluids channelled to fore-arc mantle corner: ETS and silica deposition[END_REF]. At deeper levels, water released by the dehydrating subducting plate will lower the solidus temperature, triggering the partial melt of mantle material that will migrate upward to feed the volcanic arc (S. A. [START_REF] Peacock | Fluid processes in subduction zones[END_REF].

Despite the important role played by fluids in arc magmatism and more generally in the processes occurring in subduction zones, fluid transfers and pathways in the mantle wedge and forearc crust remain largely elusive. The motivation of this study is to investigate the potential of full waveform inversion of teleseismic body waves to provide new tomographic constraints on deep subduction processes. In particular, recent methodological development has shown that accounting for the correlation between density and seismic velocities improves very significantly the reconstruction of density and V P /V S models [START_REF] Kan | A consistent multiparameter Bayesian full waveform inversion scheme for imaging heterogeneous isotropic elastic media[END_REF]. These two parameters are of special interest because they are reliable indicators of the presence of water, partial melt (e.g., [START_REF] Nakajima | Three-dimensional structure of v p , v s , and v p /v s beneath northeastern Japan: Implications for arc magmatism and fluids[END_REF][START_REF] Watanabe | Effects of water and melt on seismic velocities and their application to characterization of seismic reflectors[END_REF], or serpentinized mantle [START_REF] Christensen | Serpentinites, peridotites, and seismology[END_REF].

For the first application of our multi-parameter full waveform inversion approach, we focus on Cascadia, which stands among the most well-studied subduction zone worldwide and represents the end member of hot and young subduction, where extensive shallow dehydration of the subducting plate is expected to occur. Fluid processes are thus expected to be strongly expressed in the structural variations of seismic velocities, and in particular of the V P /V S ratio.

Along the Cascadia subduction margin, which extends from northern California to southern British Columbia, the Juan de Fuca (JdF) plate underthrusts the northern America plate with a convergence rate that ranges from 30 mm/yr at its southern end and 45 mm/yr in its northern end [START_REF] Wilson | Confidence intervals for motion and deformation of the Juan de Fuca Plate[END_REF]. In contrast to Northern Cascadia where seismicity is abundant, earthquakes in central Cascadia are scarce, so the geometry of the subducting slab [START_REF] Mccrory | Juan de Fuca slab geometry and its relation to Wadati-Benioff zone seismicity[END_REF]) is poorly constrained and tomographic studies rely primarily on seismic recordings of distant teleseismic sources. Regional travel time tomography has found that the fast anomaly associated with the subducting Juan de Fuca plate can be traced down to the transition zone except beneath central Oregon where it disappears around 200 km depth [START_REF] James | Slab fragmentation, edge flow and the origin of the Yellowstone hotspot track[END_REF][START_REF] Obrebski | Slab-plume interaction beneath the Pacific northwest[END_REF][START_REF] Obrebski | Lithosphere-asthenosphere interaction beneath the Western United States from the joint inversion of body-wave traveltimes and surface-wave phase velocities[END_REF][START_REF] Schmandt | Complex subduction and small-scale convection revealed by body-wave tomography of the Western United States upper mantle[END_REF]. The Juan de Fuca is thus fragmented (e.g., [START_REF] Obrebski | Slab-plume interaction beneath the Pacific northwest[END_REF], and the limited depth extension of the subduction beneath central Oregon which implies a small slab pull force, may explain the lack of deep seismicity in that region.

In Cascadia, as in many other subduction zones, receiver function studies have revealed the presence of a thin (less than 10 km thick) low seismic velocity and high V P /V S ratio layer [START_REF] Audet | Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing[END_REF][START_REF] Audet | Slab morphology in the cascadia fore arc and its relation to episodic tremor and slip[END_REF][START_REF] Hansen | Nature of the low velocity zone in Cascadia from receiver function waveform inversion[END_REF]; M. G. [START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF][START_REF] Bostock | The Moho in subduction zones[END_REF][START_REF] Tauzin | Multi-mode conversion imaging of the subducted Gorda and Juan de Fuca plates below the North American continent[END_REF]. Episodic tremor and slip have been well studied in Cascadia [START_REF] Rogers | Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip[END_REF]. While tremor locations are poorly constrained, especially for their depths [START_REF] Wech | Cataloging tectonic tremor energy radiation in the Cascadia subduction zone[END_REF][START_REF] Wech | Automated detection and location of Cascadia tremor[END_REF], tremors seem to be related to the overpressured oceanic crust, that is, to the thin east-dipping low-velocity layer evidenced by receiver function studies [START_REF] Audet | Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing[END_REF]. The serpentinization of the forearc mantle in Cascadia is now well documented by both geological and geophysical observations (S. M. [START_REF] Peacock | Large-scale hydration of the lithosphere above subducting slabs[END_REF][START_REF] Hyndman | Serpentinization of the forearc mantle[END_REF][START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF][START_REF] Blakely | Subduction-zone magnetic anomalies and implications for hydrated forearc mantle[END_REF]. Local earthquake tomography studies in northern Cascadia have documented very low V P /V S (lower than 1.7) in the forearc crust [START_REF] Audet | Possible control of subduction zone slow-earthquake periodicity by silica enrichment[END_REF][START_REF] Guo | Correlation of porosity variations and rheological transitions on the southern Cascadia megathrust[END_REF][START_REF] Ramachandran | The fate of fluids released from subducting slab in northern Cascadia[END_REF][START_REF] Peacock | High pore pressures and porosity at 35 km depth in the Cascadia subduction zone[END_REF][START_REF] Hyndman | Cascadia subducting plate fluids channelled to fore-arc mantle corner: ETS and silica deposition[END_REF][START_REF] Savard | Seismicity, metamorphism, and fluid evolution across the northern Cascadia fore arc[END_REF]. Such low values can only be explained by a significant quartz enrichment [START_REF] Chevrot | The Poisson ratio of the Australian crust: Geological and geophysical implications[END_REF]. This would suggest that part of the silica-saturated fluids expelled by dehydration of the subduction slab reach the crust where the temperature decrease will produce massive quartz mineralizations. This process may have important consequences for the formation and evolution of a silica-enriched continental crust.

Here, we invert the waveforms of nine teleseismic events recorded by the stations of the CASC93 experiment in order to obtain models of density, V P , V S , and V P /V S beneath central Oregon, from the surface down to 300 km depth. We first describe the data selection and preparation, and the principles of the inversion method. We then present the results of the inversion and discuss the main features of the final model. Finally, we provide interpretations of the structural anomalies in density, seismic velocities, and V P /V S in terms of the different subduction processes that occur at depth as well as comparisons with previous geophysical studies. Positions of the CASC93 broadband stations (green triangles) and physiographic provinces of Oregon. The gray triangles show the stations that were discarded owing to the absence of well-recorded teleseismic events. The blue dashed line outlines the inversion domain (Figure 5). The color circles correspond to the M ≥ 3.0 earthquakes in the USGS catalog from 1990 to 2022. The red triangles show the positions of Holocene volcanoes. The black dashed lines represent the isovalues of slab depth, taken from the Slab2 model [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF]. CR: Coastal Range, WV: Willamette Valley, WC: Western Cascades, DB: Deschutes Basin.

𝛿𝛿𝛿𝛿𝑖𝑖 = log 10 ( ∫ 𝑡𝑡 2 𝑡𝑡 1 |𝑢𝑢𝑖𝑖(𝑡𝑡)| 2 𝑑𝑑𝑡𝑡 ∫ 𝑡𝑡 2 𝑡𝑡 1 |𝑢𝑢(𝑡𝑡)| 2 𝑑𝑑𝑡𝑡 ) . ( 1
)
where u i (t) is the observed trace at station i and 𝐴𝐴 𝑢𝑢(𝑡𝑡) is the reference average wavelet. For P waves, these residuals are measured on the vertical component, whereas for SH waves, they are measured on the transverse component. When visualized in map view, the travel time and amplitude residuals allow us to quickly identify stations with abnormal time or amplitude, which are removed from the data set.

The data selection procedure is illustrated with event 9 (Figure 2). The waveforms are filtered between 10 and 25 s periods, which corresponds to the frequency band used in the first inversion. The vertical component of station A28 and A30 is flat, with very small amplitudes compared with those at the other stations. A large negative travel time anomaly is observed at station A37, presumably due to a clock problem (Figure 2a). Finally, the relative amplitude anomaly at station A16 is 0.24, which means the amplitude at this station is 1.6 times larger than at adjacent stations. We thus decided to also remove the records of this station for this particular event.

Owing to the limited duration of the acquisition, we only found 7 P and 4 SH teleseismic arrivals of sufficient quality for FWI at periods down to 4 s (Figure 3). The event information is given in Table 1.

Source Wavelet Estimation

To compute synthetic seismograms of teleseismic body waves we need the apparent incident source wavelet for each selected event. To determine these source wavelets, we first compute the Green's functions for each event that corresponds to a short duration (here 2 s) Gaussian source wavelet [START_REF] Monteiller | On the validity of the planar wave approximation to compute synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF]. These computations are performed with AxiSEM in the ak135 reference Earth model [START_REF] Kennett | Constraints on seismic velocities in the Earth from traveltimes[END_REF]. The source wavelet is then obtained by first deconvolving the Green's function from the observed vertical component at each station, using a simple water-level deconvolution in the frequency domain. The deconvolved traces are then aligned with respect to the first trace and averaged. Finally, the source wavelet is shifted so that the average time residual between synthetic and observed traces is zero. This time correction is necessary to remove the contribution of heterogeneities located outside the tomographic grid, similar to the practices in classical regional travel time tomography [START_REF] Aki | Determination of the three-dimensional seismic structures of the lithosphere[END_REF].

Design of the Tomographic Mesh

We build a spectral-element mesh that represents a chunk of the spherical Earth, centered at -122.2° longitude and 44.4° latitude (Figure 5). The size of the mesh is 6° in longitude and 2° in latitude, with hexahedral elements of about 6 km in size. Surface topography from the ETOPO1 Global Relief Model is applied on the upper edge of the spectral-element mesh. The mesh goes from the surface down to 300 km depth, the depth extension of the subducting slab in regional tomographic models beneath central Oregon (e.g., [START_REF] James | Slab fragmentation, edge flow and the origin of the Yellowstone hotspot track[END_REF][START_REF] Obrebski | Slab-plume interaction beneath the Pacific northwest[END_REF]. We use degree-7 Lagrange polynomials in the three orthogonal directions so that the geometry of each element is defined by 8 × 8 × 8 = 512 nodes. The grid is composed of 35,112 spectral elements and 12,622,584 nodes or grid points. The minimum resolved period of this mesh is about 2 s. Our goal is to invert the values of density, V P , and V S at each node of the spectral-element mesh which will also be our tomographic grid.

Definition of the Cost Function

We invert the vertical and radial components of P waveforms, and the transverse component of SH waveforms. The time windows used in the inversions start 5 s before the main phase arrival and end 30-60 s after the phase (P or SH), depending on the event. These time windows are sufficiently long to include not only the direct phases but also all the reverberations and conversions on the main lithospheric interfaces.

The waveform inversion problem consists in finding the model m that minimizes the cost function [START_REF] Tarantola | Inverse problem theory[END_REF]

𝜒𝜒(𝐦𝐦) = 1 2 (𝐠𝐠(𝐦𝐦) -𝐮𝐮obs) 𝑡𝑡 𝐂𝐂 -1 𝐃𝐃 (𝐠𝐠(𝐦𝐦) -𝐮𝐮obs) + 𝜆𝜆 2 ( 𝐦𝐦 -𝐦𝐦 prior ) 𝑡𝑡 𝐂𝐂 -1 𝐌𝐌 ( 𝐦𝐦 -𝐦𝐦 prior ) , ( 2
)
where g is the forward wave equation operator, C D the data covariance matrix, and C M the a priori model covariance matrix. The first term in Equation 2 quantifies the waveform misfits, where integration over selected time windows and summation over events and seismic phases is implicit. The second term quantifies the model deviation from the a priori model. The regularization parameter λ is introduced so that these two terms have balanced contributions to the misfit function.

The Full Waveform Inversion Method

The full waveform inversion method follows the general workflow described in [START_REF] Monteiller | Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM-DSM hybrid method[END_REF] and [START_REF] Kan | A consistent multiparameter Bayesian full waveform inversion scheme for imaging heterogeneous isotropic elastic media[END_REF], to which the reader is referred for further details. The forward wave-equation problem is solved with a hybrid numerical method [START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF][START_REF] Monteiller | On the validity of the planar wave approximation to compute synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF] that couples AxiSEM [START_REF] Nissen-Meyer | AxiSEM: Broadband 3-D seismic wavefields in axisymmetric media[END_REF] with a spectral-element method (SEM) in a regional domain [START_REF] Komatitsch | Introduction to the spectral-element method for 3-D seismic wave propagation[END_REF]. The principle is to compute and inject the incident wavefield on the edges of the SEM grid in order to limit the time-consuming 3D computations inside a small regional domain, thereby reducing the computational costs by several orders of magnitude.

Whereas it is common practice to consider a diagonal model covariance matrix, that is, to assume that model parameters are uncorrelated, [START_REF] Kan | A consistent multiparameter Bayesian full waveform inversion scheme for imaging heterogeneous isotropic elastic media[END_REF] have demonstrated that such assumption leads to inconsistent and 1.

biased model reconstructions, and that accounting for the well-documented correlations between density and V P , and between V P and V S , improved substantially the reconstructed density and velocity models. Following their approach, we consider a model covariance matrix of the form

𝐂𝐂𝐌𝐌 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 𝜎𝜎 2 𝜌𝜌 𝐂𝐂 𝑟𝑟𝜎𝜎𝜌𝜌𝜎𝜎𝑉𝑉 𝑃𝑃 𝐂𝐂 𝑟𝑟𝜎𝜎𝜌𝜌𝜎𝜎𝑉𝑉 𝑆𝑆 𝐂𝐂 𝑟𝑟𝜎𝜎𝜌𝜌𝜎𝜎𝑉𝑉 𝑃𝑃 𝐂𝐂 𝜎𝜎 2 𝑉𝑉 𝑃𝑃 𝐂𝐂 𝑟𝑟𝜎𝜎𝑉𝑉 𝑃𝑃 𝜎𝜎𝑉𝑉 𝑆𝑆 𝐂𝐂 𝑟𝑟𝜎𝜎𝜌𝜌𝜎𝜎𝑉𝑉 𝑆𝑆 𝐂𝐂 𝑟𝑟𝜎𝜎𝑉𝑉 𝑃𝑃 𝜎𝜎𝑉𝑉 𝑆𝑆 𝐂𝐂 𝜎𝜎 2 𝑉𝑉 𝑆𝑆 𝐂𝐂 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ( 3
)
where r is the correlation coefficient between each pair of parameters, C the spatial covariance kernel, σ ρ the standard deviation of density, 𝐴𝐴 𝐴𝐴𝑉𝑉 𝑃𝑃 the standard deviation of V P , and 𝐴𝐴 𝐴𝐴𝑉𝑉 𝑆𝑆 the standard deviation of V S . For the spatial covariance of a given parameter between two points of the model located at r and r′, we consider an exponential kernel of the form

𝐶𝐶 ( 𝐫𝐫, 𝐫𝐫 ′ ) = 𝑒𝑒 - √ |𝐫𝐫-𝐫𝐫 ′ | 2 𝓁𝓁 , ( 4
)
where ℓ is the correlation (or smoothing) length. Its main advantage is that it admits the simple inverse operator

𝐶𝐶 -1 ( 𝐫𝐫, 𝐫𝐫 ′ ) = 1 8𝜋𝜋𝓁𝓁 3 ( 𝐈𝐈 -𝓁𝓁 2 Δ ) 2 , ( 5
)
where Δ is the Laplacian operator [START_REF] Oliver | Calculation of the inverse of the covariance[END_REF][START_REF] Trinh | Bessel smoothing filter for spectral-element mesh[END_REF]. The values of the hyperparameters used in the inversions are given in Table 2.

The Iterative Inversion Algorithm

We search for the minimum of ( 2) with an iterative algorithm (e.g., [START_REF] Tarantola | Inverse problem theory[END_REF], in which the model m k at iteration k is perturbed by the quantity

Δ𝐦𝐦𝑘𝑘 = -𝛼𝛼𝑘𝑘𝐇𝐇 -1 𝑘𝑘 ⋅ 𝛾𝛾𝑘𝑘 ( 6
)
where α is the step length, H -1 the inverse Hessian, and γ k the gradient of cost Equation 2 at iteration k. This gradient is given by

𝛾𝛾𝑘𝑘 = 𝜕𝜕𝜕𝜕𝑘𝑘 𝜕𝜕𝐦𝐦 = 𝐆𝐆 𝑡𝑡 𝑘𝑘 𝐂𝐂 -1 𝐃𝐃 (𝐠𝐠(𝐦𝐦𝑘𝑘) -𝐮𝐮obs) + 𝜆𝜆𝐂𝐂 -1 𝐌𝐌 ( 𝐦𝐦𝑘𝑘 -𝐦𝐦prior ) , ( 7
)
with G k the derivative or Jacobian operator given by

𝐆𝐆𝑘𝑘 = 𝜕𝜕𝐠𝐠 𝜕𝜕𝐦𝐦 (𝐦𝐦𝑘𝑘). (8)
The first contribution to the gradient is computed with the adjoint method [START_REF] Tromp | Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels[END_REF]. Instead of computing the inverse Hessian, we use its finite-difference approximation given by the L-BFGS method [START_REF] Nocedal | Numerical optimization[END_REF] from the model and gradients obtained at the previous iterations. We refer the reader to [START_REF] Monteiller | Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM-DSM hybrid method[END_REF] for a more detailed description of the L-BFGS method and its application to waveform inversion.

The hierarchical inversion starts at long period with data bandpass filtered between 25 and 10 s. We then run new inversions on data filtered with a lower corner period that is progressively decreased from 10 to 8, 6, and finally 4 s, using as a starting model the model obtained at the previous longer period inversion.

Event Normalization

Our teleseismic data set contains events with different magnitudes, ranging between 6.2 and 7.5 (Table 1). We thus need to normalize the gradients of each event so that all the events contribute equally to model reconstruction. Since each event gradient is obtained by correlating the direct and adjoint wavefields [START_REF] Tromp | Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels[END_REF], we define a normalization that scales with the product of the characteristic amplitudes of the direct and adjoint wavefields. Because the adjoint wavefield is the solution of the wave equation with the adjoint sources that are

Journal of Geophysical Research: Solid Earth

KAN ET AL.

10.1029/2022JB025506

7 of 21 the time-reversed waveform residuals at each station, it scales with the root mean square waveform residual. The normalization N k for event k is thus defined by

𝑁𝑁𝑘𝑘 = 1 𝑁𝑁(𝑡𝑡2 -𝑡𝑡1) √ √ √ √ 𝑁𝑁 ∑ 𝑖𝑖=1 ∫ 𝑡𝑡 2 𝑡𝑡 1 𝑠𝑠 2 𝑖𝑖 (𝑡𝑡) 𝑑𝑑𝑡𝑡 ⋅ √ √ √ √ 𝑁𝑁 ∑ 𝑖𝑖=1 ∫ 𝑡𝑡 2 𝑡𝑡 1 [ 𝑠𝑠 𝑖𝑖(𝑡𝑡) -𝑢𝑢𝑜𝑜𝑜𝑜𝑠𝑠𝑜𝑖𝑖(𝑡𝑡) ] 2 𝑑𝑑𝑡𝑡𝑜 (9)
where s i (t) and u obs,i (t) are respectively the synthetic and observed seismograms at station i. For P waveforms, we use the vertical component to compute the normalization factors whereas for SH waveforms we use the transverse component. These normalization factors are computed at the beginning of each inversion and remain the same during iterations.

Results

Figure 6 shows the evolution of the normalized cost function during the inversions in the different frequency bands.

The misfit reduction is about 80% in the long-period inversion. Note that the misfit reductions are not directly comparable because they are normalized by the misfit of the starting model, which is not the same in the different inversions.

Tomographic Model Obtained After Hierarchical Inversion

The final density, V P , and V S models obtained after hierarchical inversion are shown in Figure 7, using both absolute and relative scales. The most salient feature is the eastward dipping low-velocity zone (LVZ) on the western side of the profile, with compressional velocities between 5.6 and 6.2 km/s and shear velocities between 3.2 and 3.5 km/s. The top of this LVZ, observed between 20 and 40 km depth, coincides with the contour of the subducted JdF slab [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF]. Since the minimum period in the inversions is 4 s, its thickness is not well resolved, but we can estimate it to be less than 10 km.

The LVZ is also characterized by a rather large V P /V S ratio ∼1.9. Even though this feature is on the margin of what is confidently resolved by our inversion, such a large V P /V S ratio is in rather good agreement with estimation coming from the modeling of receiver function which suggests even larger values, that can reach up to 2.2 [START_REF] Audet | Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing[END_REF]; S. M. [START_REF] Peacock | High pore pressures and porosity at 35 km depth in the Cascadia subduction zone[END_REF]. The LVZ has a shallower dipping angle in its eastern termination, which is consistent with the results obtained in the inversion of short-period scattered body waves [START_REF] Rondenay | Multiparameter two-dimensional inversion of scattered teleseismic body waves 3. Application to the Cascadia 1993 data set[END_REF]; M. G. [START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF]. We note that the horizontal extension of the LVZ is strongly correlated with the epicentral location of tremors given by Wech (2021) (Figure 14). From 40 to 60 km depth, seismic velocities as well as density progressively increase downdip and reach values of about 8 km/s for V P , 3.5 km/s for V S , and 3.5 g • cm -3 for density.

A similar though less pronounced velocity increase is observed in the subducting oceanic mantle, which has a clear signature beneath 50 km depth. The upper limit of this high-velocity body again closely follows the contour of the JdF slab.

To the east of the section, that is, beneath the High Cascades and Deschutes Basin, sharp contrasts in density and V P , and to a lesser extent in V S , are observed across a nearly flat Moho, located at about 35 km depth. This is in excellent agreement with the results of the inversion of short-period scattered teleseismic body waves [START_REF] Rondenay | Multiparameter two-dimensional inversion of scattered teleseismic body waves 3. Application to the Cascadia 1993 data set[END_REF]; M. G. [START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF]. In contrast, beneath the Willamette Valley and the western part of the west ern Cascades, a smooth, weak gradient in density and seismic velocities is observed at the crust-mantle boundary. This reduced impedance contrast explains the disappearance of the Moho in the migrated sections beneath the forearc [START_REF] Rondenay | Multiparameter two-dimensional inversion of scattered teleseismic body waves 3. Application to the Cascadia 1993 data set[END_REF][START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF][START_REF] Tauzin | Multi-mode conversion imaging of the subducted Gorda and Juan de Fuca plates below the North American continent[END_REF] and in the wide-angle reflection sections [START_REF] Trehu | Crustal architecture of the cascadia forearc[END_REF]. It results from a combination of high density and seismic velocities in the forearc crust and low density and seismic velocities in the forearc mantle wedge.

Active source seismic experiments have documented high V P in the forearc, reaching 6.5 km/s at 11 km depth [START_REF] Trehu | Crustal architecture of the cascadia forearc[END_REF], in rather good agreement with our V P model. These high velocities have been attributed to a block of oceanic crust, the Siletz terrane, accreted to North America ∼50 Ma ago. The results of these active source experiments also demonstrated southward crustal thickening in the Siletz terrane, from 6 km north to at least 20 km south in central Oregon [START_REF] Trehu | Crustal architecture of the cascadia forearc[END_REF]. Such a large crustal thickness could be related to an oceanic plateau formed above a hotspot near the Farallon-Kula plate boundary [START_REF] Trehu | Crustal architecture of the cascadia forearc[END_REF][START_REF] Wells | Geologic history of Siletzia, a large igneous province in the Oregon and Washington coast range: Correlation to the geomagnetic polarity time scale and implications for a long-lived Yellowstone hotspot[END_REF]. The forearc mantle wedge is characterized by a density around 3.0 g • cm -3 , V P around 7.0 km/s, and V S around 3.8 km/s, which represents a reduction of about 15% compared to the values of a typical continental lithospheric mantle. This suggests a high degree of serpentinization which can be estimated at about 30%, considering the experimental data from [START_REF] Christensen | Elasticity of ultrabasic rocks[END_REF] and [START_REF] Christensen | Serpentinites, peridotites, and seismology[END_REF]. This value is about half of the indirect estimate obtained from the analysis of scattered teleseismic waves, which is of the order of 50%-60% (M. G. [START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF]. At shallow depth (<10 km), we image three low-velocity regions in the accretionary prism, the Willamette Valley and the Deschutes Basin, which have a particularly clear signature in the V S model.

The V P /V S ratio model displays several interesting features. The western Cascades crust is characterized by a V P /V S ratio lower than 1.65. Such low values, which are well below those of typical crustal rocks (e.g., [START_REF] Brocher | Empirical relations between elastic wavespeeds and density in the Earth's crust[END_REF][START_REF] Christensen | Poisson's ratio and crustal seismology[END_REF], have been previously documented in local earthquake tomography studies of the Cascadia forearc crust [START_REF] Guo | Correlation of porosity variations and rheological transitions on the southern Cascadia megathrust[END_REF][START_REF] Ramachandran | The fate of fluids released from subducting slab in northern Cascadia[END_REF][START_REF] Savard | Seismicity, metamorphism, and fluid evolution across the northern Cascadia fore arc[END_REF].

In contrast, beneath the Willamette Valley and part of the Coastal Range, that is, in the Siletz terrane [START_REF] Wells | Geologic history of Siletzia, a large igneous province in the Oregon and Washington coast range: Correlation to the geomagnetic polarity time scale and implications for a long-lived Yellowstone hotspot[END_REF], the V P /V S ratio is comprised between 1.7 and 1.8, which suggests a more mafic composition. Near the surface and below the coastal range, V P /V S is larger than 1.9, as expected for sediments saturated in fluids in the accretionary wedge. The sediments in the Willamette Valley have also a clear signature in the V P /V S model, in good agreement with the local earthquake tomography performed by [START_REF] Dunham | Local earthquake tomography of the central Oregon forearc using a large-N, short duration, nodal array[END_REF]. Beneath the volcanic arc, a shallow prominent positive V P /V S anomaly is imaged, that could be related to a magma reservoir. Finally, a prominent high V P /V S anomaly larger than 1.9 is observed, around 75 km beneath the volcanic arc. This region is also characterized by reduced density and seismic velocities, which could suggest the presence of melt.

Waveform Fits

The P and S waveform fits for event 9 are shown in Figure 8. We compare the observed (black lines) and synthetic (red lines) seismograms computed in the initial smooth 1-D model and in the final 3-D model (Figure 7). For P waves, the fit of the vertical component before inversion is much better than the radial component. This results from the fact that the vertical component of teleseismic P waves is dominated by the direct P wave, which is highly coherent at the regional scale. It is this very property that allows us to make robust estimate of the average incoming P wavelet. The P reverberations, such as the PmP phase reflected on the Moho, generally have small amplitudes, of the order of a few percent compared to the directly transmitted P wave. The vertical component misfits thus mainly reflect the small amplitude and travel time residuals of the P wave at each station. The radial component, on the other hand, is characterized by a more evenly distributed energy in the selected time window. Secondary arrivals in the P wave coda, such as P-to-S conversions and reflections on seismic discontinuities, have amplitudes comparable to the direct P wave on the radial component. On the western part of the transect, the signature of eastward dipping interfaces, most likely related to the LVZ, is clearly observed. Because the initial model is a smooth 1-D model, the radial component synthetic seismograms computed in this model provide a poor fit of observed data. The transverse component of SH waves is also dominated by the direct SH wave, with strong later-arriving reverberations. Note that the time window selected for the inversion of SH waveforms is shorter than the one for P waveforms, in order to exclude the ScS and SKS waves that arrive ∼40 s after the S wave. After inversion, the improvement of P and SH waveforms fits is spectacular (>80%), not only for the energetic P and SH waves but also for the later arriving scattered waves.

Synthetic Resolution Test

To estimate the resolution potential of our FWI, we performed a synthetic inversion experiment on the subduction model from [START_REF] Kan | A consistent multiparameter Bayesian full waveform inversion scheme for imaging heterogeneous isotropic elastic media[END_REF]. In this study, they considered a regular 2D array of receivers with a 10 km spacing and four teleseismic sources in backazimuth of 0°, 70°, 180°, and 300°. This new synthetic inversion experiment allows studying the limitation resulting from the 1D CASC93 acquisition. The models obtained after hierarchical inversion of the synthetic seismograms computed for the P and SH teleseismic sources are shown in Figure 9 down to 4 s period.

The main features of the subduction model such as the sedimentary basin, the subducting plate, and the serpentinized mantle wedge are correctly reconstructed. However, the V P /V S model is contaminated by significant artifacts even though the main features of the input model are distinguishable. This simple inversion experiment thus suggests that the CASC93 experiment data set should indeed provide robust models of density, V P , and V S down to 150 km depth.

Discussion

The Impact of Adding SH Waveforms and Accounting for the Correlation Between Model Parameters

To investigate the impact of adding SH waveforms and of the a priori assumption that model parameters are correlated, we performed three additional hierarchical inversions:

• Inversion of only P waveforms with a diagonal model covariance matrix, that is, assuming that model parameters are uncorrelated • Inversion of both P and SH waveforms with a diagonal model covariance matrix • Inversion of only P waveforms with a complete model covariance matrix

The waveform fits for event 9 obtained after these three hierarchical inversions and the ones already presented in the previous section are shown in Figure 10. The reduction of the vertical and radial P waveform misfits for this event is similar in the four inversions. The main difference is observed for the SH waveform fits, which are greatly improved when SH waveforms are included in the inversions. Inverting only the P waveforms with a complete covariance matrix yields intermediate fits for the SH waveforms (40% misfit reduction), suggesting that accounting for the correlation between V S and density and between V S and V P partially compensates for the absence of SH waveforms for the V S model reconstruction.

Figure 11 gives the waveform misfits computed on the complete data set. We again observe similar P waveform fits in the four inversions, with a slighter larger misfit reduction on the vertical component compared to the radial component. A good SH waveform fit is only obtained when SH waves are included in the inversion.

Whereas the four inversions provide quasi-identical P waveform fits, the final V P models differ substantially (Figure 12). The models obtained with a diagonal model covariance matrix, that is, assuming that density and seismic velocities are not correlated, only show long wavelength structures such as the crust and its lateral thickness variations, as well as the subducting slab. Structural details such as the LVZ, which has a strong signature on the radial component of P waves (Figure 7c), is not retrieved at all. On the other hand, it is seen in the density models, as well as in the V S models, but with a reduced amplitude of the low-velocity anomaly. Intuitively, this can be understood by the fact that for a given frequency content, the wavelength of P waves is about twice that of S waves, and since in FWI the resolution approximately scales with the wavelength (e.g., [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF]) a sharper resolution is indeed expected in the V S model. Another reason is that the V P model is mainly reconstructed from transmitted P waves, whereas for the density and V S models the dominant contribution comes from P-to-S conversions and reflections. The former are mainly sensitive to long wavelength structures whereas the latter are mainly sensitive to sharp gradients and discontinuities. A consequence of the uneven resolution in the V P and V S models is that the V P /V S model is poorly constrained, and with strong artifacts, notably in the forearc, where extreme values lower than 1.50 are observed (Figures 12d and12h).

These results underscore the non-unique solutions of FWI and the intrinsic trade-offs between density, V P , and V S . Small waveform misfits therefore cannot be interpreted as a reliable indicator of the quality of the inversion results. Other informed criteria describing our degree of knowledge on model parameters need to be considered. Our results confirm that accounting for the well-documented correlations between V P and V S , and between density and V P (or V S ) lead to better-resolved models, especially of V P and V S models at depth. The common assumption (often made implicitly) that model parameters are not correlated, has a strong detrimental effect on the results of FWI.

The addition of SH waveforms, which were not considered in previous applications of FWI on real data (Y. [START_REF] Wang | The deep roots of the Western Pyrenees revealed by full waveform inversion of teleseismic P waves[END_REF]Beller, Monteiller, Operto, et al., 2018;[START_REF] Wang | Full-waveform inversion of high-frequency teleseismic body waves based on multiple plane-wave incidence: Methods and practical applications[END_REF][START_REF] Chevrot | Passive imaging of collisional orogens: A review of a decade of geophysical studies in the Pyrénées[END_REF], has also a beneficial though less noticeable effect. Whereas large-scale structures in the V S models shown in Figures 12k and12o are extremely similar, small-scale structures such as the LVZ, the serpentinized mantle wedge, or the low-velocity anomaly around 75 km depth beneath the volcanic arc, are more sharply defined in the inversion that includes SH waveforms.

An important limitation of FWI of teleseismic P waves comes from reverberations in sedimentary basins. In the earlier applications of FWI in the Pyrenees region, strong artifacts were indeed observed in the crust and shallow mantle beneath the Aquitaine Basin (Y. [START_REF] Wang | The deep roots of the Western Pyrenees revealed by full waveform inversion of teleseismic P waves[END_REF][START_REF] Wang | Full-waveform inversion of high-frequency teleseismic body waves based on multiple plane-wave incidence: Methods and practical applications[END_REF][START_REF] Chevrot | Passive imaging of collisional orogens: A review of a decade of geophysical studies in the Pyrénées[END_REF]. Similar artifacts are observed in the eastern part of the section, beneath the Deschutes Basin, where values of V S larger than 4.0 are observed in the lower crust (Figure 12k). Delph et al. ( 2018) also reported very high shear velocities in the crust from their joint inversion of surface wave dispersion curves and receiver functions. In contrast, values more typical of a continental lower crust (between 3.6 and 3.8 km/s) are obtained after joint inversion of P and SH waveforms (Figure 12o). Therefore, transmitted shear waves seem to provide sufficient independent information to correct the velocity model from the detrimental effects of reverberations on shallow interfaces such as the basement of sedimentary basins, which leads us to advocate the exploitation of shear waves in future applications of FWI. [START_REF] Beller | Probing depth and lateral variations of upper-mantle seismic anisotropy from full-waveform inversion of teleseismic body-waves[END_REF] came to a similar conclusion regarding the importance of jointly inverting compressional and shear waveforms for proper 3-D imaging of seismic anisotropy. We can thus foresee a growing interest in performing such joint P and S inversions in the near future.

Comparison Between FWI, Travel Time Tomography, and Scattered Waves Inversion

Figure 13 shows the comparison between our V S model and some published V S models obtained by regional travel time tomography in a vertical cross-section that follows the CASC93 transect line. The geometry of the long-wavelength fast velocity anomaly which corresponds to the JdF subducting plate is fairly consistent in all these models. However, the FWI model appears to contain a wealth of details at shallow depth (<100-150 km that are completely absent in the models derived from travel time tomography). This marked difference in resolution between FWI and travel time tomography stems from the fact that the latter only exploits the phase of the direct P or S phases. Therefore, only long-wavelength structures are correctly resolved. Moreover, since teleseismic body waves have subvertical angles of incidence, tomographic images suffer from poor sensitivity in the shallow part of the model where there is little ray crossing as well as strong vertical smearing. In contrast, FWI exploits the complete seismic wavefield and the contribution of scattered phases strongly increases the sensitivity to small-scale heterogeneities and to density. However, owing to geometrical spreading, the sensitivity of scattered waves to structural discontinuities decreases as depth increases. It is also limited by the finite duration of the temporal windows that are inverted, which can potentially suppress the contribution of waves scattered on very deep discontinuities. In the long-period FWI, the dominant period of seismic waves is longer than 10 s. This means that the wavelength of P waves is greater than 80 km and the wavelength of S waves is greater than 50 km. Therefore, the asymptotic description of the wavefield in terms of well-isolated direct or scattered P and S waves is no longer valid up to ∼100 km depth [START_REF] Favier | Near-field influences on shear wave splitting and traveltime sensitivity kernels[END_REF]. In this depth interval, the near-field and mid-field terms are not negligible. The transition to the asymptotic far-field regime at depth also corresponds to a sharp decrease in the amplitude of the gradient (see Figure 6 in [START_REF] Monteiller | Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM-DSM hybrid method[END_REF]). This sharp decrease in sensitivity where the distance Figure 4. Illustration of the source wavelet estimation procedure on event 9 (see Table 1). The vertical and transverse components are used for P and SH waveforms, respectively. (a) Observed vertical components filtered in the period band 10-25 s. (b) Time-aligned traces (gray lines) after being deconvolved from the corresponding Green's functions and the average source wavelet (red line). (c, d) Same as (a, b) but for the transverse component SH waveforms. The average time shift Δt is used to correct the timing of the source wavelet so that for each event the average travel time residual is zero. from seismic stations exceeds about one wavelength results from the fact that in the far-field regime, the main contribution comes from directly transmitted waves, with much weaker contributions from scattered waves.

Because the number of teleseismic events with a sufficiently high S/N ratio is generally small, especially when dealing with short-duration deployments such as the CASC93 Experiment, the ray coverage at depth is quite poor, in contrast to travel time tomographic studies that often rely on ∼10 5-6 travel time picks. For example, in their tomography of the mantle beneath the United States [START_REF] Schmandt | P and S wave tomography of the mantle beneath the United States[END_REF], used 516,688 P and PKP traveltimes, and 223,462 S and SKS travel times.

In Figure 14, we superpose our V S model on the results obtained by the inversion of short-period teleseismic P waves coda, which exploits the different scattering modes to reconstruct the geometry of seismic discontinuities in the underlying medium [START_REF] Rondenay | Multiparameter two-dimensional inversion of scattered teleseismic body waves 3. Application to the Cascadia 1993 data set[END_REF][START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF]. A remarkable agreement is observed with our velocity model, in particular for the geometry of the LVZ and of the Moho in the upper plate.

The absence of a clear Moho signature was attributed to the presence of higher seismic velocities in the lower forearc crust and lower velocities in the serpentinized mantle wedge (M. G. [START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF]. This interpretation is in excellent agreement with our 3D V S model. Migration of scattered waves also documented a weaker signature of the oceanic crust below ∼50 km depth [START_REF] Rondenay | Multiparameter two-dimensional inversion of scattered teleseismic body waves 3. Application to the Cascadia 1993 data set[END_REF][START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF][START_REF] Tauzin | Multi-mode conversion imaging of the subducted Gorda and Juan de Fuca plates below the North American continent[END_REF], which is also in good agreement with our FWI results. These studies also detected a dipping interface characterized by a downward velocity reduction between 50 and 100 km depth beneath the oceanic Moho. This interface has been attributed to the upper boundary of a layer of frozen anisotropy resulting from the alignment of olivine crystals near the ridge [START_REF] Tauzin | Multi-mode conversion imaging of the subducted Gorda and Juan de Fuca plates below the North American continent[END_REF]. We find no evidence for a such velocity reduction in our model.

The Geometry of the Subducting Juan de Fuca Plate

The Slab2 subduction zone model [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF] is built from a compilation of seismicity catalog and tomographic models. In Cascadia, due to the paucity of the earthquake catalog, the JdF plate is primarily constrained by published travel time tomography models such as those shown in Figure 13. Although these different models differ in detail, they show consistent long-wavelength structures with a fast anomaly down to about 200 km depth approximately following the slab contour from the Slab2 model. They also show a pronounced low-velocity anomaly beneath the backarc between 100 and 300 km depth. Our FWI model is in good agreement with these first-order features. The increase of slab dip around 50 km depth corresponds to the completion of eclogitization reactions in the oceanic crust. These dehydration reactions are associated with an increase in the density of the oceanic crust that makes the plate negatively buoyant. Eclogitization of the oceanic crust could thus explain the apparent change in slab dip. In contrast, the results of scattered wave inversion [START_REF] Rondenay | Multiparameter two-dimensional inversion of scattered teleseismic body waves 3. Application to the Cascadia 1993 data set[END_REF][START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF][START_REF] Tauzin | Multi-mode conversion imaging of the subducted Gorda and Juan de Fuca plates below the North American continent[END_REF], based on the negative polarity interface imaged between 50 and 100 km, favor a smaller dip of about 20°, close to that of the LVZ. With such a geometry, the low seismic velocities observed around 100 km depth in our model (Figure 14) could be related to the partial melting of the subducting oceanic crust. Numerical modeling experiments suggest that partial melting of a subducting oceanic crust is indeed possible in subduction zones where very young oceanic lithosphere is being subducted, such as in the Cascade Range (S. M. [START_REF] Peacock | Partial melting of subducting oceanic crust[END_REF]. Although this hypothesis cannot be excluded, we tend to favor the first hypothesis which is more consistent with tomographic models of the Cascade Range. Regarding the disappearance of the fast anomaly below 200 km depth, since our model has a very limited spatial resolution in this depth range, our study does not shed new light on this long-standing controversy.

Fluid Migrations, Melting, and Metamorphic Reactions in the Central Cascadia Subduction Revealed by FWI

Our tomographic image (Figure 7) reveals the 3-D structures of the Cascadia subduction in central Oregon with a wealth of details that provide crucial insight into fluid migration and related metamorphic mineral reactions in the subduction factory. Our interpretations of the model are summarized in Figure 13.

The nature of the ubiquitous low-velocity layer beneath the Cascadia forearc, which is a prominent feature of our tomographic model, has been discussed by M. [START_REF] Bostock | The Moho in subduction zones[END_REF]. It is thought to be composed of hydrated sediments and pillow basalts. Dehydration reactions release fluids that remain trapped between a low permeability seal above the subduction interface and an underlying layer of low-porosity gabbroic rocks. Tremors in central Cascadia are weaker and less frequent than in northern or southern Cascadia [START_REF] Wech | Cataloging tectonic tremor energy radiation in the Cascadia subduction zone[END_REF]. Episodic tremor and slip (ETS) [START_REF] Rogers | Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip[END_REF] locations from the Pacific Northwest Seismic Network (https://pnsn.org/ tremor) [START_REF] Wech | Cataloging tectonic tremor energy radiation in the Cascadia subduction zone[END_REF] are closely related to the LVZ [START_REF] Gosselin | Seismic evidence for megathrust fault-valve behavior during episodic tremor and slip[END_REF], which suggests a hydrologic control of slow slip events and tremors. The depth of tremors is not well constrained, especially in central Oregon where the station distribution is rather sparse. They probably occur at the subduction interface (black stars in Figure 15), as do the low-frequency events to which they are strongly related [START_REF] Thomas | Identifying low-frequency earthquakes in central Cascadia using cross-station correlation[END_REF]. However, a depth distribution of tremors over a depth range exceeding 40 km has been documented beneath southern Vancouver Island [START_REF] Kao | A wide depth distribution of seismic tremors along the northern cascadia margin[END_REF]. Relocating tremors in finely resolved structural models such as the one obtained in this study could potentially improve our understanding of the fluid processes that produce these enigmatic events.

In our model, we observe a strong contrast between the backarc mantle and the forearc mantle, where density and seismic velocities are reduced by about 15% compared to the backarc mantle. The limit between these two domains is quite sharp and approximately follows the geometry of isotherm 600° in the thermal structure calculated for a hot continental subduction [START_REF] Hyndman | Serpentinization of the forearc mantle[END_REF]; Figure 4a). The transition from a low-permeability layer interface to a high-permeability layer interface corresponds to the onset of eclogitization of the oceanic crust [START_REF] Audet | Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing[END_REF]. The eclogitization domain is clearly observed in our model, between the LVZ and the eclogitic crust, where the density and seismic velocities are close to those of the surrounding mantle. The width of this eclogitization domain, which is about 40 km, allows us to estimate the kinetics of the reaction. The velocity of the JdF subducting plate is about 4 cm/year [START_REF] Demets | Current plate motions[END_REF] and thus a complete eclogitization of the oceanic crust takes ∼1 Myr. This is an important conclusion since eclogitization reactions (e-h) Inverted models using the CASC93 stations and the selected events in Table 1.

(i-l) Inverted models using the source-receiver configuration considered in [START_REF] Kan | A consistent multiparameter Bayesian full waveform inversion scheme for imaging heterogeneous isotropic elastic media[END_REF], with a regular 2-D array and four teleseismic sources in backazimuth of 0°, 70°, 180°, and 300°. Each column shows models of absolute density, V p , V s , and V p /V s from the first row to the fourth row.

Figure 10. Waveform fits for event 9 obtained after the four hierarchical inversions. The solid lines show the observed seismograms and the red lines show the synthetic seismograms. The blue lines show the SH synthetic seismograms that were obtained in the two inversions in which only P waveforms were inverted. The cost reduction in percentage for each component after hierarchical inversion is specified at the upper-right of each panel.

control the buoyancy of subducting plates. Such a rapid reaction rate suggests that the subduction of an oceanic plateau is unlikely to support flat subduction, as suggested by geodynamic modeling studies [START_REF] Gerya | Dynamic effects of aseismic ridge subduction: Numerical modelling[END_REF][START_REF] Van Hunen | On the role of subducting oceanic plateaus in the development of shallow flat subduction[END_REF].

The transformation of hydrated metabasalts to eclogite releases a large amount of water into the overlying mantle wedge. This water will hydrate the olivine minerals in the peridotites, producing serpentine minerals [START_REF] Hyndman | Serpentinization of the forearc mantle[END_REF]. Serpentine and hydrous minerals significantly alter the physical and mechanical properties of the forearc mantle, with density and seismic velocities reduced in proportion to the degree of serpentinization [START_REF] Christensen | Serpentinites, peridotites, and seismology[END_REF]. Back-arc mantle is too hot for hydrous minerals to be stable [START_REF] Hyndman | Serpentinization of the forearc mantle[END_REF].

Serpentinization thus leads to a weak and buoyant forearc mantle wedge, isolated from the subduction-related circulation flow. The eastward limit of the serpentinized mantle wedge is quite sharp, and approximately corresponds to the 600°C isotherm in [START_REF] Hyndman | Serpentinization of the forearc mantle[END_REF]. As already noted by M. G. [START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF], this also suggests a significant temperature contrast between the forearc and the backarc mantle, which could explain the sharp eastward increase of heat flow measured at the surface [START_REF] Blackwell | Heat flow, arc volcanism, and subduction in northern Oregon[END_REF]. Serpentinization reactions also produce magnetite. Since the temperature in the forearc mantle wedge is expected to be cooler than 580°C [START_REF] Hyndman | Serpentinization of the forearc mantle[END_REF], the Curie temperature of magnetite, a magnetic anomaly is thus expected along the Cascadia forearc. Such a large anomaly is indeed observed, associated with a gravity low [START_REF] Blakely | Subduction-zone magnetic anomalies and implications for hydrated forearc mantle[END_REF] and thus the hypothesis of a highly hydrated mantle wedge is also supported by potential field anomalies.

Fluid pathways in the mantle wedge are still unclear, but tomographic studies suggest that fluids do not migrate vertically (e.g., [START_REF] Halpaap | Seismicity, deformation, and metamorphism in the Western Hellenic subduction zone: New constraints from tomography[END_REF][START_REF] Schurr | Complex patterns of fluid and melt transport in the central Andean subduction zone revealed by attenuation tomography[END_REF]. Instead, the fluids follow oblique diverging paths oriented toward the volcanic arc and/or upward the subduction interface toward the accretionary wedge, separated by a zone in the mantle wedge with low-fluid content. These intriguing observations were recently explained by numerical modeling of water migration in the mantle wedge with a two-phase flow approach (H. [START_REF] Wang | Water migration in the subduction mantle wedge: A two-phase flow approach[END_REF]. The water migration models in this study are characterized by the development of two water pathways with a low-porosity divide in between, which result from the bending of the overlying lithosphere that produces a zone of overpressure (which can reach up to 50% of the lithostatic pressure) and low-porosity in the mantle wedge. If we relate the decrease of density and seismic velocities observed in the serpentinized mantle wedge to the amount of mantle hydration, our FWI model clearly suggests two diverging fluid pathways with a core that has higher density and seismic velocities, located beneath the eastern margin of the Willamette Valley. Magnetotelluric studies imaged conductive features that approximately correspond to these two hypothesized fluid pathways [START_REF] Evans | Electrical structure of the central Cascadia subduction zone: The EMSLAB Lincoln Line revisited[END_REF][START_REF] Wannamaker | Segmentation of plate coupling, fate of subduction fluids, and modes of arc magmatism in Cascadia, inferred from magnetotelluric resistivity[END_REF]. A first prominent conductive anomaly (resistivity <1 Ω•m) which diverges eastward, can be traced from the top of the slab to the Moho beneath the western Cascades, from where it extends horizontally further to the east. Interestingly, [START_REF] Vidale | Deep long-period earthquakes west of the volcanic arc in Oregon: Evidence of serpentine dehydration in the fore-arc mantle wedge[END_REF] reported deep long-period earthquakes (red circles in Figure 15) located ∼40 km west of the volcanic arc, near the Moho, that is, coinciding with the horizontal part of the first conductivity anomaly. Deep long-period earthquakes are commonly observed beneath volcanoes, associated with the movement of fluids [START_REF] Chouet | A multi-decadal view of seismic methods for detecting precursors of magma movement and eruption[END_REF]. The location beneath the western Cascades near the Moho thus gives further support to the hypothesis of fluids rising from the slab, percolating through the serpentinized mantle wedge, and reaching the base of the forearc crust. The second conductive anomaly follows the plate interface updip and emerges offshore, at the front of the accretionary wedge.

Around 100 km depth, fluids released by dehydration of the slab reduce the solidus temperature and produce partial melting. The presence of melt will lower the shear velocity and, to a lower extent, the compressional velocity, and thus a positive V P /V S anomaly is a good indicator of partial melt. Our model indeed reveals an anomaly that possesses all the attributes of partial melt, beneath the volcanic arc at around 70 km depth, which could indeed be the deep source of the magma plumbing system. Magnetotelluric data are equivocal on whether there is a conductor that could be related to partial melt at this depth, because of the mask produced by the strong shallow conductive anomaly beneath the arc [START_REF] Evans | Electrical structure of the central Cascadia subduction zone: The EMSLAB Lincoln Line revisited[END_REF].

The forearc crust is characterized by an anomalously low V P /V S ratio, which can only by explained by large-scale silica enrichment [START_REF] Brocher | Empirical relations between elastic wavespeeds and density in the Earth's crust[END_REF][START_REF] Chevrot | The Poisson ratio of the Australian crust: Geological and geophysical implications[END_REF][START_REF] Christensen | Poisson's ratio and crustal seismology[END_REF]. This enrichment requires massive transport of fluids saturated in silica that migrate upward through the forearc crust and precipitate when the solubility of silica drops because of the temperature decrease [START_REF] Manning | The solubility of quartz in H 2 O in the lower crust and upper mantle[END_REF]. Massive and pervasive mineralized quartz veins are indeed a well-documented ubiquitous feature of forearc crusts (e.g., [START_REF] Breeding | Slab-derived fluids and quartz-vein formation in an accretionary prism, Otago Schist, New Zealand[END_REF][START_REF] Halpaap | Seismicity, deformation, and metamorphism in the Western Hellenic subduction zone: New constraints from tomography[END_REF][START_REF] Ramachandran | The fate of fluids released from subducting slab in northern Cascadia[END_REF], which may have contributed to the silica enrichment of continents. These quartz veins often host gold deposits which result from extensive leaching of the crust by hot fluids [START_REF] Kerrich | Nature's gold factory[END_REF].

Conclusions

We have imaged the deep architecture of the subduction beneath central Cascadia, from the surface to 300 km depth, by inverting complete P and SH teleseismic waveforms. Compared to classical tomographic studies, the teleseismic FWI approach simultaneously constrains density, V P , V S , and the V P /V S ratio with fine spatial resolution. The final tomographic models provide new insight into the dehydration of the subducted Juan de Fuca slab and the associated eclogitization metamorphic reactions. They also reveal the fluid pathways through the forearc mantle as well as a partial melt zone around 75 km depth beneath the volcanic arc.

The main limitations of this study stem from the linear geometry of the acquisition, and from the short duration of the CASC93 experiment. We thus expect that more detailed and robust tomographic images could be obtained in northern Cascadia, where the stations are both more numerous and better distributed. Another important perspective will be to extend full waveform inversion to general anisotropic media [START_REF] Beller | Probing depth and lateral variations of upper-mantle seismic anisotropy from full-waveform inversion of teleseismic body-waves[END_REF] to constrain the 3D distribution of seismic anisotropy in the Cascadia subduction zone. [START_REF] Dunham | Local earthquake tomography of the central Oregon forearc using a large-N, short duration, nodal array[END_REF] and deep long-period events [START_REF] Vidale | Deep long-period earthquakes west of the volcanic arc in Oregon: Evidence of serpentine dehydration in the fore-arc mantle wedge[END_REF], respectively. The thick black line draws the subduction interface that corresponds to the downgoing Juan de Fuca plate. Straight thick blue arrows indicate the path of ascending fluids originating from the slab that react with the olivine minerals of the mantle wedge to produce serpentinized mantle with reduced density and seismic velocities. Wiggled thin blue arrows indicate the hypothesized paths followed by fluids in the forearc crust. The Moho beneath the forearc and the volcanic arc is defined from the 3.8 km/s isocontour in the V S model. The 600°C isotherm taken from [START_REF] Hyndman | Serpentinization of the forearc mantle[END_REF] is indicated by a thin black line.

Figure 1 .

 1 Figure1. Positions of the CASC93 broadband stations (green triangles) and physiographic provinces of Oregon. The gray triangles show the stations that were discarded owing to the absence of well-recorded teleseismic events. The blue dashed line outlines the inversion domain (Figure5). The color circles correspond to the M ≥ 3.0 earthquakes in the USGS catalog from 1990 to 2022. The red triangles show the positions of Holocene volcanoes. The black dashed lines represent the isovalues of slab depth, taken from the Slab2 model[START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF]. CR: Coastal Range, WV: Willamette Valley, WC: Western Cascades, DB: Deschutes Basin.

Figure 2 .

 2 Figure 2. Example of data selection procedure for the vertical component of event 9 in Table 1. (a) Travel time residuals. Negative residuals represent phases arriving earlier than the average trace. (b) Amplitude residuals. (c) Correlation coefficient of each station and the average trace. (d) Observed vertical data filtered between 10 and 25 s period.

Figure 3 .

 3 Figure 3. Teleseismic events for (a) P waves (b) and SH used in this study. Detailed event informations are given inTable 1.

Figure 5 .

 5 Figure 5. Spectral-element mesh. The number of elements along the x, y, and z directions are 56, 19, and 33, respectively. The yellow triangles show the positions of the stations from the CASC93 experiment that were used in this study. The black contour shows the coastline. The background colors in the grid show the initial V S model, a smoothed version of the ak135 reference Earth model. Note the velocity model near the boundary of the mesh is specified for the layered ak135 model, to be consistent with the one used to compute the external wavefield.

Figure 6 .

 6 Figure 6. Evolution of cost function reductions for the four hierarchical inversions with progressively reduced lower corner period of 10, 8, 6, and 4 s. The misfits are normalized to 100% when we start the new inversion in a shorter period.

Figure 7 .

 7 Figure 7. Results of the hierarchical full waveform inversion. Left: models of absolute (a) density, (b) V p , (c) V s , and (d) V p /V s . Right: perturbations of (e) density, (f) V p , (g) V s , and (h) V p /V s with respect to the initial 1-D smooth ak135 model. The black line shows the slab contour of the JdF plate taken from Hayes et al. (2018), and the red triangle marks the position of Mount Jefferey in the volcanic arc. Physiographic provinces are the same as specified in Figure 1.

Figure 8 .

 8 Figure8. Waveform fits for event 9 at filtered between 4 and 25 s period (see Table1). Observed (black lines) and synthetic (red lines) vertical component P waveforms computed in the initial (a) and final (b) models. Observed (black lines) and synthetic (red lines) radial component P waveforms computed in the initial (c) and final (d) models. Observed (black lines) and synthetic (red lines) transverse component SH waveforms computed in the initial (e) and final (f) models. The black ticks on each trace indicate the beginning and end of the time windows used in the inversion. The reductions of cost functions for the vertical, radial, and transverse components for this event are 85%, 81%, and 89%, respectively. The secondary arrivals observed ∼40 s after the S wave in (e,f) are the SKS and ScS phases.

Figure 9 .

 9 Figure 9. Reconstructions of the synthetic subduction model. (a-d) True model. (e-h) Inverted models using the CASC93 stations and the selected events in Table1. (i-l) Inverted models using the source-receiver configuration considered in[START_REF] Kan | A consistent multiparameter Bayesian full waveform inversion scheme for imaging heterogeneous isotropic elastic media[END_REF], with a regular 2-D array and four teleseismic sources in backazimuth of 0°, 70°, 180°, and 300°. Each column shows models of absolute density, V p , V s , and V p /V s from the first row to the fourth row.

Figure 11 .

 11 Figure 11. Residual misfits for the different types of inversions, defined as the ratio between the misfit in the final and the misfit in the initial smooth 1-D model (in percents).

Figure 12 .

 12 Figure12. Results of the hierarchical FWI with different inversion strategies. Each column shows models of absolute density, V p , V s , and V p /V s from the first row to the fourth row. (a-d) Inversion with P events and a diagonal covariance matrix. (e-h) Inversion with P and SH events and a diagonal covariance matrix. (i-l) Inversion with P events and a non-diagonal covariance matrix. (m-p) Same as Figures7a-7d. Inversion with P and SH events and a non-diagonal covariance matrix.

Figure 13 .

 13 Figure 13. Vertical cross-sections showing V S models along the CASC93 experiment transect (latitude of 44.44°N). The dashed line shows the slab contour of the JdF plate (Hayes et al., 2018). (a) FWI V S model from this study. The rectangular box represents the region shown in Figure 7. (b) V S model from DNA10 (Obrebski et al., 2011). (c) V S model from NWUS11-S (James et al., 2011). (d) V S model from DNA13 (Porritt et al., 2014). (e) V S model from US-SL-2014 (Schmandt & Lin, 2014). (f) V S model from Casc19-S (Hawley & Allen, 2019).

Figure 14 .

 14 Figure 14. Comparison between our final V S model and the main interfaces reconstructed by inversion of short period scattered teleseismic body waves (M. G.[START_REF] Bostock | An inverted continental Moho and serpentinization of the forearc mantle[END_REF]. The tremor intensity represents the number of tremors in every 2 km bin with a horizontal search radius of 20 km. Note the coincidence between the horizontal extent of the low-velocity zone and the occurrence of tremors.

Figure 15 .

 15 Figure 15. Interpretation of the V S model shown in Figure 7, with (vertical exaggerated) topography plotted on top. The green circles show the positions of the stations of the CASC93 temporary experiment. The gray and red circles indicate the locations of crust earthquakes (Dunham & Kiser, 2020) and deep long-period events (Vidale et al., 2014), respectively. The thick black line draws the subduction interface that corresponds to the downgoing Juan de Fuca plate. Straight thick blue arrows indicate the path of ascending fluids originating from the slab that react with the olivine minerals of the mantle wedge to produce serpentinized mantle with reduced density and seismic velocities. Wiggled thin blue arrows indicate the hypothesized paths followed by fluids in the forearc crust. The Moho beneath the forearc and the volcanic arc is defined from the 3.8 km/s isocontour in the V S model. The 600°C isotherm taken from Hyndman and Peacock (2003) is indicated by a thin black line.

Table 1

 1 List of Teleseismic Events Selected for Waveform Inversion

	Event	Phases	Origin time	Location (lat, lon)	Depth (km)	Dist. (°)	Baz (°)	M w
	1	P	1993/05/16 21:44:50.0	15.34 S, 173.35 W	33	75.4	231	6.6
	2	SH	1993/07/06 02:53:03.9	24.56 S, 111.86 W	10	69.3	170	6.3
	3	P	1993/07/11 13:36:19.2	25.35S, 70.18 W	33	84.2	134	6.6
	4	P	1993/07/12 13:17:11.9	42.84 N, 139.25 E	17	66.7	308	7.3
	5	P + SH	1993/08/07 17:53:27.0	23.87 S, 179.82 E	555	86.1	231	6.7
	6	P	1993/08/08 08:34:25.0	12.96 N, 144.78 E	61	83.1	282	7.5
	7	SH	1993/08/11 14:17:41.7	13.13 N, 145.62 E	56	82.3	281	6.2
	8	P	1993/09/10 19:12:54.8	14.73 S, 92.68 W	34	38.9	130	7.2
	9	P + SH	1993/10/11 15:54:22.4	32.00 N, 137.85 E	365	74.6	300	6.8

Table 2

 2 Hyperparameters Used in the Inversions

	Periods 𝐴𝐴	𝐴𝐴𝑉𝑉 𝑃𝑃 (km/s) 𝐴𝐴	𝐴𝐴𝑉𝑉 𝑆𝑆 (km/s)	σ ρ (g/cm 3 )	r	ℓ (km)
	25-10, 8, 6 s	0.65	0.37	0.27	0.97	5
	25-4 s	0.65	0.37	0.27	0.80	5

21699356, 2023, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JB025506 by Portail BibCNRS INSU, Wiley Online Library on [11/04/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Acknowledgments

We thank Roy Hyndman and an anonymous reviewer for their insightful and constructive reviews. V. Monteiller was supported by the European Union's Horizon 2020 research and innovation programme under the ChEESE project, grant agreement No. 823844. This work was performed using HPC resources from GENCI-IDRIS (Grant 2021-101085) and CALMIP (Grant P20012). The authors also acknowledge PRACE for awarding us access to Marconi100 at CINECA, Italy.

Data Availability Statement

The teleseismic waveforms are publicly available from Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC). The FDSN code of the CASC93 experiment is XZ (1993XZ ( -1994)). The FWI code is available upon request to the corresponding author, and it will be released to the public soon.


Djamal Habet
email: djamal.habet@ema.fr

Chu Min Li

Laure Devendeville

Michel Vasquez
email: michel.vasquez@ema.fr

A Hybrid Approach for SAT

Exploiting variable dependencies has been shown very useful in local search algorithms for SAT. In this paper, we extend the use of such dependencies by hybridizing a local search algorithm, Walksat, and the DPLL procedure, Satz. At each node reached in the DPLL search tree to a fixed depth, we construct the literal implication graph. Its strongly connected components are viewed as equivalency classes. Each one is substituted by a unique representative literal to reduce the constructed graph and the input formula. Finally, the implication dependencies are closed under transitivity. The resulted implications and equivalencies are exploited by Walksat at each node of the DPLL tree. Our approach is motivated by the power of the branching rule used in Satz that may provide a valid path to a solution, and generate more implications at deep nodes. Experimental results confirm the efficiency of our approach.

Introduction

Consider a propositional formula F in Conjunctive Normal Form (CNF) on a set of boolean variables {x 1 , x 2 , . . . , x n }, the satisfiability problem (SAT) consists in testing wether all clauses in F can be satisfied by some consistent assignment of truth values to variables.

SAT is the first known [START_REF] Cook | The Complexity of Theorem Proving Procedures[END_REF] and one of the most well-studied NP-complete problem. It has many applications like graph coloring, circuit designing or planning, since such problems can be encoded into CNF formulas in a natural way.

Stochastic Local Search (SLS) approaches for SAT became prominent, when independently Selman, Levesque, and Mitchell [START_REF] Selman | A New Method for Solving Hard Satisfiability Problems[END_REF] as well as Gu [START_REF] Gu | Efficient Local Search for Very Large-Scale Satisfiability problems[END_REF] introduced algorithms based on stochastic local hill-climbing. They are considered as the most powerful incomplete methods for solving large and hard SAT instances. However, such algorithms can get stuck in the local minima of the search space, and do not integrate the structural relations between variables in their resolution. On the other hand, the complete methods, based on the Davis-Putnam-Logemann-Loveland procedure (DPLL) [START_REF] Davis | A Machine Program for Theorem Proving[END_REF], depend on the choice of the variable to branch on. One of the best recent implementations of DPLL procedure,

Walksat

Originally introduced in [START_REF] Selman | Noise Strategies for Improving Local Search[END_REF], Walksat performs a greedy local search for a satisfying assignment of a set of propositional clauses in SAT format. The procedure starts with a randomly generated truth assignment. It then changes (flips) the assignment of a variable chosen under the heuristic described in algorithm 2 below. Flips are repeated until either a satisfying assignment is found or a preset maximum number of flips, Max-Flips, is reached. This process is repeated up to a maximum number of Max-Tries times. Compared to its predecessors, like GSAT, Walksat differs in one important aspect. In fact, while GSAT architecture is characterized by a static neighborhood relation between assignments with Hamming distance one, Walksat's one is based on a dynamically determined subset of the GSAT neighborhood relation. Effectively, the experimental results show that Walksat outperforms the existing SLS algorithms proposed before, and it is proved Probabilistically Approximately Complete (PAC property) with a noise setting wp > 0 [START_REF] Hoos | Local search algorithms for SAT: An empirical evaluation[END_REF].

Satz

Despite its simplicity and seniority, the Davis-Putnam-Logemann-Loveland procedure (DPLL) remains one of the best complete procedures for SAT. It essentially constructs a binary search tree and its nodes are results of recur-sive calls. While a solution is not found, all leaves represent a dead-end where a contradiction (empty clause) is found. DPLL procedure performance is closely related to the selection of the branching variable. In fact, this selection affects the search tree size, and consequently the required time to solve F . A popular and a cheap branching heuristic is the MOM1 heuristic [START_REF] Freeman | Hard Random 3-SAT Problems and the Davis-Putnam Procedure[END_REF], which picks the variable that occurs the most often in the minimal size clauses. However, work realized in Posit [START_REF] Freeman | Improvements to Propostional Satisfiability Search Algorithms[END_REF][START_REF] Freeman | Hard Random 3-SAT Problems and the Davis-Putnam Procedure[END_REF], Tableau [START_REF] Crawford | Experimental Results on the Crosover Point in Random 3-SAT[END_REF], and Satx [START_REF] Li | Exploiting Yet More the Power of Unit Clause Propagation to solve 3-SAT Problem[END_REF] have suggested to integrate unit propagations to the heuristic H. It results in an UP 2 heuristic which examines the variable x by respectively adding the unit clauses x and x into F , and independently making two unit propagations. The real effect of those propagations is then used to weigh x. However, since examining variables by two unit propagations is time consuming, it is necessary to reduce the number of variables examined by the UP heuristic. Taking into account the number of binary occurrences of variables, Satz gets the best restrictions on the number of examined variables. In that way, combining the MOM and UP heuristics, Satz reduces the size of the search tree by detecting failed literals as early as possible.

Hybrid Approach

Variable Dependencies

SAT encodings of structured problems, such as planning and diagnosis, often contain large numbers of variables whose values are constrained to be a simple boolean function of other variables. These variables are then dependent. Variables whose values cannot be easily determined to be a simple function of other variables are independent. For a given SAT problem, there may be many different ways to classify the variables as dependent. In this work, we use the dependency definition given in [START_REF] Lang | Complexity Results for Independence and Definability in Propositional Logic[END_REF] as follows:

Definition 1. Let Σ be a set of clauses, V the related variables, C finite con- junctions of literals such that V C ⊆ V ,

and a variable y such that Σ |= (y ⇒ C).

Then we say that the variables of V C depend on the variable y. Roughly speaking, if y is instantiated then all variables of V C are instantiated too. The implications are naturally constructed by unit propagations, performed by Satz when looking for the variable to branch on, and are represented by a directed graph, where nodes are literals and edges relate two dependent literals. This construction is done in linear time.

Favorably, two literals mutually implied are equivalent. An equivalency is a stronger dependency than a simple implication. A set of equivalent literals, constituting a class, is a strongly connected component of the implication graph, where a representative literal, chosen randomly, substitutes the other element of its class. Consequently, the input formula is reduced and the implication graph is also reduced and becomes acyclic.

Dependencies Consistency

In the preprocessing of a SAT formula or when a branching occurs, Satz fixes some variables and their states become passive. To maintain the implications coherent, this state must be propagated through the implication graph for the satisfied literals, and through its transposed graph for the falsified literals. Those propagations are performed by a Depth-First-Search procedure in a linear time, which is also used to construct the equivalency classes restricted to the active variables. On the other hand, the implication relation is transitive. Once again, by the Depth-First-Search procedure, the implication closure is constructed from the reduced implication graph. At this stage, the implications are fully enumerated and are consistent.

Proposed Algorithm: WalkSatz

After applying the enumerated processes to each sub-formula obtained in a current node of the Satz tree, Walksat is then applied by considering the equivalency classes rather than original variables. The class to flip is selected as described in the section 2. The implications are to be integrated advantageously in the Walksat resolution. In fact, three cases are distinguished when Walksat chooses a variable to flip: a zero-damage-flip, a minimal-damage-flip, and a random walk. It is then possible to merge the implication constraint in one, two, or all the three cases, i.e., in any case, if a class is flipped then the implied ones are flipped too. We have tested the eight possible combinations under a large variety of problems, and the best integration was observed in the case of a minimal-damage-flip. Such result can be interpreted by: 1. The zero-damage-flip is the best case, but flipping the implied classes may increase the violated clause number. 2. The random walk is incompatible with the deterministic behavior of implications, and integrating the implication relation to all levels, without any improved mechanism, the search may be trapped easily in the local minima.

We should remark that a large reduction of an instance size may cause, in the undesirable cases, a cancellation of successive actions (flips) applied to a same class. To avoid this, a tabu list [START_REF] Glover | Tabu Search[END_REF] of a tenure fixed empirically to 1, is used as follows: each flipped class, as well as the related ones by the implication relation (if flipped), are forbidden to any change during one iteration. The full steps of the hybrid approach are resumed in the developed algorithm, WalkSatz, which is incomplete because of the restriction on the depth of the Satz tree. The hybrid approach is the result of combining a systematic method and a local search one. Consequently, it is not so evident to give comparative criteria between such families. However, in order to evaluate the performances of WalkSatz, we compare it with Walksat and Satz on a broad range of benchmarks Despite the importance of the noise parameter, no optimization is made on. However, its value is fixed to an identical value, 0.5, to both algorithms. The depth of the Satz tree is limited to 3, the root being at level 0, at most fifteen nodes are then developed, and the total number of the authorized flips is identical for WalkSatz and Walksat. The columns T%, # flips, and sec. correspond respectively to the success rate, average flips, and average run time in second. # flips and sec. are calculated for the successful executions, and the running time for Satz is limited to 7200 seconds. All the used instances are satisfiable.

Latin Square

Given a set S, a Latin square indexed by S is an |S| × |S| array such that each row and each column of the array are a permutation of the elements in S. |S| is called the order of the Latin square. These instances have been contributed by H. Zhang [START_REF] Zhang | Implementing the davis-putnam method[END_REF].

DIMACS Benchmarks

"aim" instances: proposed by Iwama and al. [START_REF] Asahiro | Random Generation of Test Instances with Controlled Attributes[END_REF], the instances are all generated with a particular random 3-SAT instance generator.

"ssa" instances: contributed by A.V. Gelder and al., the instances correspond to single-stuck-at-faults problem in circuit analysis. The used instances are selected formulas from those generated by N emesis [START_REF] Ferguson | Test Pattern Generation for Realistic Bridging Faults in CMOS ICS[END_REF][START_REF] Larrabee | Test Pattern Generation Using Boolean Satisfiability[END_REF], which is a testpattern generation program. "par8" instances: contributed by J. Crawford and suggested by M. Kearns, these instances are propositional versions of the parity learning problem.

Superscalar Processor Verification

Defined by M.N. Velev [START_REF] Velev | Superscalar processor verification using efficient reductions of the logic of equality with uninterpreted functions to propositional logic[END_REF], these instances result from the verification of exceptions, multicycle functional units, and branching prediction in superscalar microprocessors, DLX.

Beijing-Challenge Benchmarks

This set comprises the instances which have been proposed for the International Competition on SAT Testing in Beijing, 1996, including planning and scheduling problems [START_REF] Kautz | Encoding Plans in Propositional Logic[END_REF][START_REF] Kautz | Pushing the Envelope: Planning, Propositional Logic, and Stochastic Search[END_REF].

Kautz & Selman Planning Problems

Proposed by H. Kautz and B. Selman [START_REF] Kautz | Encoding Plans in Propositional Logic[END_REF][START_REF] Kautz | Pushing the Envelope: Planning, Propositional Logic, and Stochastic Search[END_REF], these instances correspond to the planning logistic problem.

Experimental Result Discussion

Globally and when Compared to Walksat, WalkSatz largely reduces the required flips number to get a solution, and has a good behavior when solving hard instances. These statements may be justified as follows:

1. The equivalency classes reduce the problem size and eliminate the cyclic implications. On the other side, the implication transitive closure permit to enumerate all the dependent literals, and so, exploit fully their advantage. 2. The implication relation can be viewed as a tabu which restricts the visited neighborhood by the local search. Furthermore, checking the implication consistency allows to verify implicit constraints, then reduce the number of flips required to reach a solution. 3. The branching rule used by Satz gives two sub-formulas completely disjoint.

So, the space search explored by Walksat is diversified.

Excepting "qg1-08", "qg5-11", and "4blocks" instances, the high rate success of WalkSatz shows the robustness of this hybrid algorithm. Furthermore, whenever Satz works well, WalkSatz takes advantage in robustness, and when Satz works less well, for example for the "dlx2 cc" instances, there is a small improvement in robustness at the expense of computation time. On random SAT problems containing few variable dependencies, Walksatz and Walksat essentially have the same behavior. The experimental results are not presented here. As said in the section 5, it is less evident to make an adjusted comparison, and the used parameters are retained to compare the behavior of WalkSatz and WalkSat with a same number of authorized flips. However, other parameters are to underline: in the above experiments, if WalkSatz fails in a given node of the Satz tree, no restart is done on (only one try is authorized), but other experiments show that WalkSatz works better with more tries: a restart is an useful aspect in a local search. The depth of the Satz tree is also prominent, and its optimized value can be found for each problem class. However, a tardeoff should be respected: it is look not so interesting to develop many nodes, in the Satz tree, to let the local search reaching a solution. At the end, the tabu tenure and list, used by WalkSatz, are to be well handled: because of the instance reduction, these tabu parameters can be a function of the number/size of the equivalency classes.

The main weakness of WalkSatz is the computation time requirements. It is mainly caused by the repeated Depth-First-Search procedures at every node of the Satz search tree. In fact, the version elaborated in this paper is not optimized. This lack is easily improved by an incremental construction of the implications and equivalencies at each branching and backtracking done by Satz. Despite this drawback, the experimental results prove the efficiency of our approach.

Conclusion

In this paper, we presented a hybrid approach between two efficient algorithms Walksat and Satz. Our work is motivated by the performances of such algorithms. In fact, by its powerful branching rule, Satz provides an optimal search tree, and by its variable selection heuristic, Walksat outperforms the precedent SLS algorithms. On the other hand, SAT formulas contain implicit variable dependencies that are so useful in a SAT resolution, especially in a local search process. To well exploit those elements, WalkSatz was elaborated. It is the result of numerous experiments in the aim to get a good hybridization that increases the contribution of implication and equivalency relationships in the resolution mechanism. The obtained results support our approach.

As perspective works, more powerful variants of Walksat, such as Novelty and R-Novelty [START_REF] Mcallester | Evidence for Invariants in Local Search[END_REF], are also indicated for such approach. The variable dependencies, as well as tabu, noise setting, and other WalkSatz parameters, are to be further studied to increase their effects. Finally, to improve the computation time, these dependencies must be constructed incrementally.

Example 1 .

 1 Consider the set of clauses Σ ={¬l ∨ a, ¬l ∨ ¬b, b ∨ ¬c}. If l is fixed then the literals a, ¬b, ¬c are fixed too. So we have the dependencies l → {a, ¬b, ¬c}. Such dependencies are implications.

Algorithm 1 :

 1 WalksatInput: SAT-formula F, Max-Tries,Max-Flips Output: A satisfying truth assignment T of F, if found

	begin
	for try=1 to Max-Tries do
	T ←randomly generated truth assignment;
	for flip=1 to Max-Flips do
	if T satisfies F then return T ;
	c←randomly selected clause violated under T ;
	v←Heuristic(F,c);
	T ←T with v flipped;
	return "Solution not found";
	end;
	Algorithm 2: Heuristic

Input: SAT-formula F, violated clause c Output: Selected variable to flip, v begin for each variable u appearing in c do Calculate score(u) equal to the violated clause number in F if u is flipped; if there are variables with null score then v←randomly select one (zero-damage-flip); else switch a probability value (noise setting) do case wp : v←variable with minimal score (minimal-damage-flip); otherwise v ←randomly select a variable from c (random walk); return v ; end;

  Construct for F, the implication graph I and its transposed I t ; Propagate variables state under I and I t ; Construct the equivalency classes and reduce F, I, and I t ; Construct the implication closure for the reduced graphs I and I t ; Apply Walksat to F reduced to its equivalency classes, taking into account the reduced graphs I and I t when minimizing the total number of violated clauses (minimal-damage-flip);

	if the last step returns "Solution T found" then
	return "Solution T , by Walksat";
	return "Solution not found";
	end;

Algorithm 5: WalkSatz Input: SAT-formula F Output: A satisfying truth assignment T of F, if found begin for each node reached by Satz, down to a fixed depth do if F is empty then return "Solution T found, by Satz";

  3 . WalkSatz is coded in Linux/C++ and compiled with g++ compiler. Comparative parameters used in our experiments 4 are:

		WalkSatz Walksat
	Max-Tries	1	100
	Max-Flips	10 5	15 × 10 5
	Number of runs	100	1
	Noise setting	0.5	0.5

Table 1 .

 1 Experimental results for the "qg" instances

	Latin square	WalkSatz		Walksat		Satz
	Instance # Vars. # Cls. T% # flips.	sec. T% # flips. sec.	sec.
	qg1-08	512 148957 6 1307271 1499.446 0	-	-	125.220
	qg2-07	343	68083 100 39719	0.901 25 693660 694.622 24.340
	qg3-08	512 148957 100 238622 1.505	9 628604 148.794 0.210
	qg4-09	512	10469 83 701706 9.169	1 328792 3226.650 0.750
	qg5-11	1331 64054 8 277167 370.769 0	-	-	4.830
	qg6-09	729	21844 100 368	0.664	0	-	-	0.830
	qg7-09	729	22060 100	81	0.824	0	-	-	0.970

Table 2 .

 2 Experimental results for the "aim" instances

	AIM		WalkSatz	Walksat		Satz
	Instance	# Vars. # Cls. T% # flips. sec. T% # flips. sec.	sec.
	aim-50-1-6-yes1-1	50	300 100	1	0.001 6 920684 44.422 0.020
	aim-50-1-6-yes1-2	50	300 100	1	< 0.001 1 64016 148.600 0.030
	aim-50-1-6-yes1-3	50	300 100	1	< 0.001 34 685077 6.378 0.010
	aim-50-1-6-yes1-4	50	300 100	1	< 0.001 0	-	-	0.010
	aim-50-2-0-yes1-1	50	100 100	1	< 0.

Table 3 .

 3 Experimental results for the "ssa" instances

	SSA	WalkSatz	Walksat	Satz
	Instance # Vars. # Cls. T% # flips. sec. T% # flips. sec. sec.
	ssa7552-160 1391 3126 100 385 2.270 100 30875 0.061 0.080
	ssa7552-159 1363 3032 100 344 2.265 100 22006 0.433 0.080
	ssa7552-158 1363 3034 100 255 2.218 100 26738 0.052 0.070
	ssa7552-038 1501 3575 100 991 2.581 100 80673 0.167 0.120

Table 4 .

 4 Experimental results for the "par8" instances

		PAR		WalkSatz	Walksat	Satz
	Instance # Vars. # Cls. T% # flips. sec. T% # flips. sec. sec.
	par8-1-c	64	254 100 12062 0.024 100 10628 0.025 0.02
	par8-1	350	1149 100 68175 0.180 4 977367 79.366 0.01
	par8-2-c	68	270 100	1	0.003 100 15504 0.036 0.02
	par8-2	350	1157 100 120942 0.356 2 483886 160.791 0.02
	par8-3-c	75	298 100 35293 0.069 100 37213 0.089 0.02
	par8-3	350	1171 100 96110 0.272 3 538312 107.811 0.02
	par8-4-c	67	266 100	1	0.000 100 41423 0.098 0.02
	par8-4	350	1155 100 22966 0.681 0	-	-	0.02
	par8-5-c	75	298 100 12843 0.027 100 34534 0.084 0.02
	par8-5	350	1171 100 77604 0.284 1 304169 330.066 0.03

Table 5 .

 5 Experimental results for the "dlx cc" instances

	DLX2 CC	WalkSatz	Walksat	Satz
	Instance	# Vars. # Cls. T% # flips. sec. T% # flips. sec.	sec.
	dlx2 cc bug01 1515 12808 100 989391 49.607 79 603523 24.857 > 7200
	dlx2 cc bug02 1515 12808 100 858251 43.818 80 517174 21.958 > 7200
	dlx2 cc bug03 1515 12808 100 1013661 51.013 71 723533 33.193 > 7200
	dlx2 cc bug14 1516 12811 100 458607 26.524 98 381390 10.122 > 7200
	dlx2 cc bug16 1516 12811 100 461816 25.555 92 417319 13.432 > 7200
	dlx2 cc bug33 1516 12798 100 636107 23.721 99 265186 6.820 > 7200
	dlx2 cc bug34 1516 12718 100 965240 43.048 98 369863 9.282 > 7200
	dlx2 cc bug38 1515 12783 100 652728 35.267 63 516764 1.404 > 7200
	dlx2 cc bug39 1482 12111 93 1141143 49.067 66 640518 33.861 > 7200
	dlx2 cc bug40 1520 12811 100 629230 33.210 98 358702 9.688 > 7200

Table 6 .

 6 Experimental results for the Beijing instances

	Beijing	WalkSatz	Walksat	Satz
	Instance	#Vars. # Cls. T% # flips. sec. T% # flips.	sec.	sec.
	3blocks		370 13732 100 26768 0.222 100 49960	0.626	0.060
	4blocks		900 59285 34 519608 102.429 0	-	-	0.060
	4blocksb		540 34199 100 264532 3.655 19 694517 210.913 0.650
	e0ddr2-10-by-5-1 19500 108887 100 246496 1067.484 0	-	-	137.820
	e0ddr2-10-by-5-4 19500 104527 98 463371 3116.833 18 1202849 132.465 185.900
	enddr2-10-by-5-1 20700 111567 97 292415 2352.466 0	-	-	112.960
	enddr2-10-by-5-8 21000 113729 100 261973 1463.593 2 1329010 1448.225 105.230
	ewddr2-10-by-5-1 21800 118607 97 252841 1119.552 0	-	-	401.460
	ewddr2-10-by-5-8 22500 123329 100 205542 846.741 1 1456023 3136.267 120.640

Table 7 .

 7 Experimental results for the Kautz & Selman planning problem

	Planning	WalkSatz	Walksat		Satz
	Instance # Vars. # Cls. T% # flips. sec. T% # flips. sec.	sec.
	logistics.a 828	6718 100 21508 0.088 100 133275 0.469 2.750
	logistics.b 843	7301 100 36783 0.142 100 249287 0.887 0.090
	logistics.c 1141 10719 100 77653 0.334 87 490317 2.896 0.420
	logistics.d 4713 21991 100 286041 155.530 74 624500 6.990 507.490
	bw large.a 459	4675 100 5146 0.199 100 14359 0.057 0.070
	bw large.b 1087 13772 100 96838 1.899 86 564590 6.249 0.240
	bw large.c 3016 50157 100 198672 22.738 0	-	-	2.740

Maximum Occurrences of clauses of Minimum size

http://www.satlib.org

All experiments are on a Duron 800 Mhz machine with 256 MB of RAM.

This work is partially supported by French CNRS under grant number SUB/2001/0111/DR16