Foivos Karakostas
email: foivos.karakostas@ingv.it

Nicholas Schmerr

Ross Maguire

Quancheng Huang

Doyeon Kim

Vedran Lekic

Ludovic Margerin

Ceri Nunn

Sabrina Menina

Taichi Kawamura

Philippe Lognonné

Domenico Giardini

Bruce Banerdt

Scattering Attenuation of the Martian Interior through Coda Wave Analysis

 for these families of events and their expected interaction with a shallow scattering layer. The shapes of the observed data envelopes do not show a consistent pattern with event distance, suggesting that the diffusivity and scattering layer thickness is non-uniform in the vicinity of InSight at Mars. Given the consistency in the envelope shapes between HF and VF events across epicentral distances and the tradeoffs between the parameters that control scattering, the dimensions of the scattering layer remain unconstrained but require that scattering strength decreases with depth and that the rate of decay in scattering strength is fastest near the surface. This is generally consistent with the processes that would form scattering structures in planetary lithospheres.

INTRODUCTION

Scattering of seismic waves from random heterogeneities is a well-studied phenomenon [START_REF] Aki | Scattering and attenuation of shear waves in the lithosphere[END_REF]; [START_REF] Ishimaru | Wave propagation and scattering in random media[END_REF]; [START_REF] Wu | Attenuation of short period seismic waves due to scattering[END_REF] and others) that depends strongly on the relative length scale of the heterogeneity and wavelength of seismic waves. For elastic waves, random heterogeneities are defined by contrasts between materials of differing seismic wave speeds, that are related to changes in shear rigidity, bulk modulus, and density of materials. Such changes are common in geologically complex materials, and expected where there are non-uniform variations in materials with depth and location. A common simplifying assumption is to use linear scaling between wavespeed and density when approximating these variations [START_REF] Sato | Unified approach to amplitude attenuation and coda excitation in the randomly inhomogeneous lithosphere[END_REF]. A further simplifying assumption in (semi)analytic approaches is that scatterers are isotropic, although full waveform methods enable the treatment of anisotropic scatterers [START_REF] Cormier | Anisotropy of heterogeneity scale lengths in the lower mantle from pkikp precursors[END_REF]. The strength of the perturbation (i.e., random heterogeneity) is thus defined as the relative change in the seismic wave speed that occurs across a discontinuous boundary in the material, usually represented with a percentile of the relevant parameter. The other aspect of heterogeneity is size; scales of random perturbations within a medium are typically characterized using an autocorrelation function, where the correlation distance is an approximation for the size of the heterogeneities within the medium (see a review of the topic by [START_REF] Shearer | Seismic scattering in the deep earth[END_REF]). In materials where the sizes of heterogeneities are large compared to the seismic wavelength, weak forward scattering dominates. If the heterogeneities are considerably larger than the seismic wavelength, then scattering effects become negligible.

Likewise, if heterogeneities are considerably smaller than the seismic wavelength, then scattering effects disappear and the medium behaves as a homogeneous solid. Scattering effects are strongest when the sizes of the scatterers and seismic wavelength are similar and there is a large seismic velocity contrast between neighboring heterogeneities [START_REF] Aki | Quantitative seismology[END_REF]. Thus, by examining the contribution of scattering effects across a range of wavelengths (or equivalently frequencies), it becomes possible to constrain the strength and size distribution of scatterers within the medium through which the waves propagate. In this study we look into a range of events with frequencies ranging from below 1 Hz to above 9 Hz. In the discussion section, we analyze the frequency-dependence of the inferred scattering properties.

The analysis of the coda of seismic body waves (e.g., P-and S-waves) at different frequencies provides valuable insight into the material properties of the medium through which they travel. This is because the rate of decay of coda energy and the signature of energy loss are related to both the intrinsic attenuation structure (1/Q i ) and seismic scattering (1/Q s ) in the subsurface. [START_REF] Wesley | Diffusion of seismic energy in the near range[END_REF] and [START_REF] Aki | Origin of coda waves: source, attenuation, and scattering effects[END_REF] used single scattering and diffusion theory to describe how Q can be inferred from S-wave coda analysis. The multi-scattering case and a diffusion equation was used by [START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda q[END_REF] in order to develop the a radiative transfer equation to study the coda waves in an inhomogeneous layered medium. Similar approaches that used the diffusion equation in order to take into account the leakage from a diffusive layer to an underlying elastic half-space were developed by [START_REF] Margerin | Residence time of diffuse waves in the crust as a physical interpretation of coda q: application to seismograms recorded in mexico[END_REF] and [START_REF] Wegler | Analysis of multiple scattering at vesuvius volcano, italy, using data of the tomoves active seismic experiment[END_REF].

In most of the Earth's interior, seismic waves experience relatively weak scattering, allowing for the direct observation of individual seismic body waves (like P-and S-waves) and surface waves. Where scattering is present, it manifests as later-arriving codas that directly follow the main arrivals and decay with time. Scattering is typically associated with small-scale compositional heterogeneities in the crust [START_REF] Revenaugh | Geologic applications of seismic scattering[END_REF], lithosphere [START_REF] Kennett | Multiscale seismic heterogeneity in the continental lithosphere[END_REF], mantle [START_REF] Mancinelli | Constraints on the heterogeneity spectrum of earth's upper mantle[END_REF], core-mantle boundary region [START_REF] Kim | Sequencing seismograms: A panoptic view of scattering in the core-mantle boundary region[END_REF][START_REF] Ma | Small-scale scattering heterogeneities in the lowermost mantle from a global analysis of pkp precursors[END_REF], and even inner core [START_REF] Leyton | Using pkikp coda to determine inner core structure: 1. synthesis of coda envelopes using single-scattering theories[END_REF] of Earth. From a wide range of studies, small-scale heterogeneities appear to be omnipresent throughout the Earth, with the exception of the well-mixed outer core, and are typically assigned to compositional variations, as thermal anomalies would not be long-lived in the Earth due to thermal diffusion [START_REF] Shearer | Seismic scattering in the deep earth[END_REF]. Although some environments like fault zones and volcanic edifices [START_REF] Prudencio | Spatial distribution of intrinsic and scattering seismic attenuation in active volcanic islands -I: model and the case of Tenerife Island[END_REF] show multiple scattering, the majority of the Earth is typically characterized by relatively weak scattering that is concentrated in the near surface, and allows for observations of direct seismic arrivals with weak codas.

In contrast, seismic waves propagating within the Moon are dominated by scattering effects (see [START_REF] Nunn | Lunar Seismology: A Data and Instrumentation Review[END_REF] for a review). A significant fraction of the seismic energy produced by moonquakes undergoing intense scattering in the near surface to the point that body waves cannot be readily identified on seismograms, and surface waves are non-existent. The lunar crust has undergone billions of years of impact gardening [START_REF] Cintala | Impact-induced thermal effects in the lunar and mercurian regoliths[END_REF]) that has produced an upper surface layer of regolith, with a seismic P-wave velocity of 100 -300 m/s [START_REF] Kovach | Apollo 14 and 16 active seismic experiments[END_REF], and underlain by a megaregolith consisting of fractured and cracked materials that may extend to 30 km depth [START_REF] Lognonné | A new seismic model of the moon: implications for structure, thermal evolution and formation of the moon[END_REF]. This impact-modified layer and a lack of intragranular fluids create a strongly-scattering environment for seismic waves. [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF] inferred the properties of a scattering layer by analytically relating them to the energy envelopes of natural and artificial impact events. They identified the density of the scatterers in the lunar shallow structure, as well as an attenuation factor of Q = 5000 for events with dominant frequencies 0.5 to 1 Hz. Since these foundational studies, the thickness of the lunar scattering layer still remains under debate, with recent estimates extending the scattering layer to 100 km or more, with a Qs ≤ 10 (Blanchette- [START_REF] Blanchette-Guertin | Investigation of scattering in lunar seismic coda[END_REF][START_REF] Blanchette-Guertin | Effects of lateral variations in megaregolith thickness on predicted lunar seismic signals[END_REF][START_REF] Garcia | Lunar seismology: An update on interior structure models[END_REF][START_REF] Gillet | Scattering attenuation profile of the Moon: Implications for shallow moonquakes and the structure of the megaregolith[END_REF]. As strong scattering in the lunar crust and megaregolith likely results from nearly ubiquitous fractures and cracks that extend into lunar rocks, as well as from the thick blanket of ejecta deposits and impact melt that persists in the shallower crust, the properties and thickness of the layer are expected to be highly variable across the surface of the Moon [START_REF] Nakamura | Shallow lunar structure determined from the passive seismic experiment[END_REF]. Although the upper layer of scattering produces seismic codas that last well over an hour and obfuscate the detection of body waves, signal processing and polarization filtering have enabled successful detections of a lunar core [START_REF] Garcia | Very preliminary reference Moon model[END_REF][START_REF] Weber | Seismic detection of the lunar core[END_REF] as well as other internal interfaces (see [START_REF] Lognonné | 10.03 -planetary seismology[END_REF] for a review).

The landing of InSight in late 2018 [START_REF] Banerdt | Initial results from the insight mission on mars[END_REF] now presents the opportunity to study the nature of seismic waves propagating within Mars. While Mars has extensive surface impact cratering, it also possesses an atmosphere and evidence of resurfacing through erosion and deposition of sediments by liquid water and lava flows [START_REF] Carr | Mars: surface and interior[END_REF]. Therefore, it would be expected that the planet might lie somewhere between the weak scattering regime present on Earth and the highly diffusive wave propagation found on the Moon.

InSight is the first ever seismometer to operate directly on the surface of Mars, with the Viking-1 and 2 landers having also brought seismometers that were both lander mounted [START_REF] Anderson | The viking seismic experiment[END_REF]. Unfortunately Viking-1's seismometer failed to uncage, and only a single seismic event was identified on the Viking-2 seismometer [START_REF] Lazarewicz | The viking seismometry[END_REF]. InSight has successfully detected over 500 marsquakes during its two first years of operation on the surface (InSight Marsquake Service, 2021), enabling studies of scattering and seismic attenuation in the Martian interior. Based on the methodology of [START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda q[END_REF], a radiative transfer model was used by [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF] to constrain, for the first time, the attenuation and scattering structure of the Martian crust. [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF] reported different diffusivities in the Martian upper crust depending on the frequency content of the examined events. Receiver function analyses of InSight data has revealed crustal layering, including a shallow, 10-km thick, low velocity layer below the InSight lander [START_REF] Knapmeyer-Endrun | Thickness and structure of the martian crust from insight seismic data[END_REF][START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF].

Additional marsquakes, available in the 5th version of the Mars Seismic Catalog (InSight Marsquake Service, 2021), offer the opportunity to examine the scattering properties of Mars more thoroughly and constrain the characteristics of the scattering attenuation in the upper crust. Therefore, in this study, we use this extended marsquake dataset and systematically explore a series of parameters that control the seismic scattering on Mars. Additionally, we further examine and analyze the frequency-dependence of diffusivity reported by [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF]. Our objective is to provide further insights and explore the limits of scattering analysis given the available marsquake dataset and the context of a single seismometer on a planet.

Martian seismic events are classified by their spectral properties [START_REF] Giardini | The seismicity of mars[END_REF] into different event types that occur across a range of dominant frequencies and distances from the InSight lander. According to [START_REF] Clinton | The marsquake catalogue from insight, sols 0-478[END_REF], when the dominant frequency range of the event is below the 2.4 Hz resonance [START_REF] Ceylan | Companion guide to the marsquake catalog from insight, sols 0-478: Data content and non-seismic events[END_REF]Kim et al., 2021a), it is considered a Low Frequency event, whereas the Broadband events include seismic energy excitation below and above 2.4 Hz. High Frequency and Very High Frequency events are dominantly above the threshold defined by the 2.4 Hz resonance.

These are the basic features of the main categories of the events that are examined in this study:

• Low Frequency (LF) events show a very rapid coda decay (less than 1 minute). The seismic catalog (InSight Marsquake Service, 2021) locates thes events to epicentral distances greater than 30 • .

• Broadband (BB) events are located at relatively large distances (more than 30 • and their coda decay is around 1.5 times longer than that of Low Frequency events. Previous studies [START_REF] Giardini | The seismicity of mars[END_REF]) located these events in the Cerberus Fossae region.

• High Frequency (HF) events are located in the 20 • -30 • distance range. However, the distance estimates depend substantially on the choice of velocity model. HF decay times are longer than those of families where lower frequencies (below 1 Hz) are dominant.

• Very High Frequency (VF) events are further classified in two families according to their epicentral distance: one, very near family closer than 20 • ; and a second, more distant family, further than 30 • . The coda decay duration of VF events is relatively long, comparable to that of HF events, and does not appear to correlate with the epicentral distance.

Teleseismic events with identifiable P and S waves are characterized by dominant frequencies f < 1 Hz, and the body waves show some degree of coda following the initial arrivals. The lack of strong scattering and detectable surface waves in these events is interpreted to imply that they must have a deeper hypocenter and therefore the recorded waves were generated in a medium where scattering is weak. Local events are associated with higher frequency content, typically having their energy as f > 2.4 Hz, and are characterized by strong codas following the P and S waves.

There are some events across all the different types that are located at an epicentral distance of 30 • . Among them, the lower frequency events (LF and BB) exhibit a different spectral character than the higher frequency ones (HF and VF) which is interpreted to result from shallow sources and wave interaction with scatterers in the low velocity layer found in the uppermost 10 km of the crust.

van Driel et al. (2021) conducted an analysis of the High Frequency seismic events detected by InSight. They modeled the spectral envelopes of High Frequency events with the Spectral Element Method (SEM) which could replicate the behavior of the HF event energy envelopes by placing a layer of strong scattering in the uppermost portion of their crustal model. They concluded that to explain the observations of high frequency seismic energy and its propagation over significant time and distance, the shallow Martian crust had to possess high Q and some degree of scattering. However, they did not attempt to directly constrain the scattering structure with their models.

In this study, we investigate the observation that Martian events dominated by lower frequency energy appear to have shorter coda decays than those with higher frequency content. It is also observed that the envelopes of the BB events for frequencies lower than 1 Hz appear to have longer coda decays than the respective of LF events. We hypothesize that the frequency signature of these events and the coda decay are associated with their ray paths through the Martian interior. [START_REF] Giardini | The seismicity of mars[END_REF] suggested the LF events are associated with rays that crossed the mantle depths and a mantle low velocity zone (described by [START_REF] Khan | Upper mantle structure of mars from insight seismic data[END_REF]), whereas the HF (and VF) event rays are primarily trapped in the diffusive and lower velocity part of the crust (see an analysis by [START_REF] Knapmeyer-Endrun | Thickness and structure of the martian crust from insight seismic data[END_REF]). This hypothesis could explain the longer coda decays for the HF and VF events. We systematically model the codas observed on the vertical component for an expanded dataset of both lower (LF, BB) and higher frequency (HF, VF) types of events, in order to quantify the scattering properties of the Martian interior. Our approach investigates the characteristics of the S-wave coda decay for these marsquake event types, and examines the tradeoffs between the layer diffusivity, thickness, velocity ratio, and background intrinsic Q of the Martian crust. We investigate how these constraints vary across event types to expand upon results identified in three previous studies of the Martian seismic attenuation [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF][START_REF] Menina | Energy envelope and attenuation characteristics of High Frequency (HF) and Very High Frequency (VF) Martian events[END_REF][START_REF] Van Driel | High frequency seismic events on mars observed by insight[END_REF], and determine the shallow diffusive structure that is consistent across event types.

DATA The Martian Seismic Events

We analyze waveforms of 21 marsquakes recorded on the vertical component of the Very Broad Band (VBB) seismometer (InSight Mars SEIS Data Service., 2019; InSight SEIS Science Team, 2019; [START_REF] Lognonné | Seis: Insight's seismic experiment for internal structure of mars[END_REF]. We select seismic events from the Martian Seismic Catalog (InSight Marsquake Service, 2021), which classifies seismic events based on their quality and frequency content. We use 19 quality B events, which are defined by the identification of either multiple clear phases but no polarization (identifiable distance but no location) or polarization but no clear phase picks (identifiable azimuth but no distance), and 2 quality A events, which are defined by both clear phase picks and polarization and therefore their distance and azimuth are identifiable [START_REF] Clinton | The marsquake catalogue from insight, sols 0-478[END_REF].

We select InSight raw data (InSight Mars SEIS Data Service., 2019; InSight SEIS Science Team, 2019) for a time window that starts 30 minutes before and ends 90 minutes after the indicated time of each seismic event in the Events Catalog InSight Marsquake Service (2021). In each seismogram we remove the instrument response through deconvolution, and then rotate the data from the modified Galperin arrangement [START_REF] Lognonné | Seis: Insight's seismic experiment for internal structure of mars[END_REF] to vertical (Z) and horizontal (North, East) components. We then taper the data using a window size of 5% of the length of the seismogram with a Tukey window.

We select the events on the basis of the clarity of the S-wave coda decay, particularly the absence of any glitches (see [START_REF] Scholz | Detection, analysis, and removal of glitches from insight's seismic data from mars[END_REF]) that would affect the examined signal and therefore contaminate our analysis. Because the seismic data recorded by SEIS include a number of peculiar signals arising from coupling between different InSight sensors and spacecraft components (Kim et al., 2021a), we manually examine each event waveform to ensure that the effect of those signal irregularity was minimal for our analysis.

The 5th version of the Seismic Catalog (InSight Marsquake Service, 2021) contains 1 LF and 1 BB, Quality A events, and 6 LF, 2 BB, 32 HF and 10 VF, Quality B events. We exclude from our analysis waveforms with glitches present or events with no discernible S-wave arrival, yielding a final dataset comprised of 21 marsquakes that occurred between the Sols 128 and 500 of the InSight operations on Mars. Our dataset includes events in all four marsquake families; 5 Low Frequency (LF), 3 Broadband (BB), 8 High Frequency (HF), and 5 Very High Frequency (VF) events, as shown in Table 1.

Spectral envelopes selection

We calculate spectral envelopes by first computing spectrograms of the vertical velocity timeseries for each seismic event using a window length of 50 s with 90% overlap.The time window for each spectrogram starts 30 minutes before and ends 90 minutes after the official time of the event listed in the Seismic Catalog. The event appears as a distinct region of higher amplitudes. We use a visual interactive tool in order to manually choose the desired lower and upper frequency of the energy envelope and the start and end of the event timeseries. The selected frequency range for each event can be found in Table 1 and the time window of the selected waveforms is shown in Figure 1. The spectral amplitude is then summed over a desired frequency range at each point in time to create the smoothed event envelope, with a sampling rate of 5 Hz.

The 1 Hz tick noise, which is discussed in the study of [START_REF] Ceylan | Companion guide to the marsquake catalog from insight, sols 0-478: Data content and non-seismic events[END_REF], affected the selection of the examined envelopes. The amplitude of this periodic signal is high enough to affect the quality of the data (see an analysis by Kim et al. (2021a)) and therefore required us to select only the portion of BB-type events with energy below this threshold. In the frequency range below 1 Hz we observe that the coda decay times of the BB events do not differ significantly from the LF events.

Examples of energy envelopes from each frequency dependent event type are shown in Figure 1. In the timeseries of the event, we select only the part of the S-wave arrival and coda decay for our analysis. We manually pick the S-wave arrival and we use the envelope until the end of the selected time window in our analysis, in order to include the whole S-wave coda decay. We note that the selected time window almost certainly includes Sg and other crustal seismic phases; however, we demonstrate later in the Methodology section that inclusion of these phases does not affect the results of our investigation.

Figure 2 presents the complete dataset that is used in this study, and identifies where several glitches that are excluded from our analysis exist. When InSight seismograms are plotted as time series, it is difficult to immediately identify glitches. However, when we analyze the signal with our visualization tool that shows the spectrogram of each event, these glitches are readily identified by a characteristic band of energy with a broadband signature that extends from the highest to lowest frequencies, in the range of f < 0.1 Hz. The low frequency energy is diagnostic of a glitch, as none of the cataloged marsquakes possess substantial energy at these low frequencies.

METHODOLOGY

We use the InSight SEIS-VBB seismogram dataset above to investigate the characteristics of scattering attenuation in the Martian crust and upper-most mantle. To study scattering, we define a structural model consisting of a diffusive layer overlying an elastic halfspace. We assume the presence of a shallow diffusive structure similar to [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF], and focus upon the strong scattering found in the near surface of planetary bodies. The underlying elastic layer in our investigation is assumed to be a half-space, which means that we do not infer the depth of the Moho or other underlying seismic interfaces [START_REF] Kim | Improving constraints on planetary interiors with pps receiver functions[END_REF][START_REF] Knapmeyer-Endrun | Thickness and structure of the martian crust from insight seismic data[END_REF]. Furthermore, we show that the characteristics of the elastic half-space do not affect our results: the only parameter of the elastic half-space that enters our analysis is its average seismic velocity, which is only used to define velocity ratio between it and the overlying, scattering layer.

The diffusive layer is assumed to be an isotropic homogeneous layer, with a unique average diffusivity and average seismic velocity. The current data from Mars are not extensive enough to enable the mapping of lateral and azimuthal variations in these parameters. As is explained in the following section, the effective thickness of the diffusive layer seen by the seismic waves depends upon the epicentral distance and the velocity ratio with the underlying elastic layer.

The average seismic ray approach

Figure 3 shows the geometry of the two layers used in our model. The diffusive layer is assumed to contain a homogeneous distribution of scatterers. Seismic waves are generated at the source and travel in all directions. The wave propagation direction from source to receiver is thus represented by a seismic ray. The seismic rays, which are shown in white, change directions when they interact with a scatterer, as the waves are either reflected, refracted, or absorbed. The reflection of the waves by the scatterers is discussed here, whereas the effect of absorption and therefore the loss of energy in the diffusive layer is summarized through a single quality factor, Q. [START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda q[END_REF] examined in detail the transition from the elastic to diffusive regime by modeling coda waves using a solution for the radiative transfer equation [START_REF] Papanicolaou | Transport equations for the stokes parameters from maxwell's equations in a random medium[END_REF]. In this study, we use the approach of an average seismic ray, from the source to the station, in order to represent the average path of the seismic rays in the diffusive and elastic layer. The direction of this average ray is not affected by the presence of the scatterers and it follows the path of a seismic wave propagating in two elastic layers with different velocities. In other words, the diffusivity in the top layer is translated into a lower apparent seismic velocity. The range of the seismic ray is defined by the velocity ratio between the two layers, which is the ratio of the apparent velocity in the top diffusive layer and the real seismic velocity in the underlying elastic layer. The seismic velocities are related to the speed of S-waves in each medium as our investigation is focused on S-wave related attenuation.

The geometry of the seismic ray

The geometry of the average seismic ray is shown in Figure 4. The diffusive layer with an apparent S-wave velocity V d overlies the elastic half-space, which has an S-wave velocity V e . In reality, V d is not constant in the crust and it changes with depth. The reader should note that the thickness of the diffusive layer, h, is exaggerated in the Figure . In addition, it should be noted that the only discontinuity shown in this schematic representation is the one between the aforementioned two layers, whereas other discontinuities of the planetary internal structure (e.g. Moho, Core-Mantle Boundary) are omitted. Therefore, the elastic layer is a space extended towards the center of the planet, but its properties are only relevant in the depth range traveled by the examined seismic ray.

In the same Figure 4, the planet's radius is noted by R, and the epicentral distance in units of length by ∆ and radians by ε. The range of the seismic ray corresponds to the projection on the surface of the planet of the ray path in the diffusive layer, and it is noted with r in length units.

On the left side, where the ray is refracted into the elastic layer, this range is divided in two equal parts, one on the source side and one on the station. On the right side, where the ray represents the reflection, the range, r, coincides with the epicentral distance, ∆. The range is expressed in radians by θ r , which is divided in two equal angles on the left (refraction) part and coincides with the epicentral distance on the right (reflection). In the refracted ray case (on the left) the incident angle is noted with i and the angle of refraction with ψ.

We use this simplified geometry to compute the range of the seismic ray in the diffusive layer, depending on the thickness of the layer and the velocity ratio V d /V e . The reflection case on the right is part of the refraction case on the left, which the refraction involving additional wave propagation in the elastic layer. The relationship among the layer thickness to the range and velocity ratio is given in Equation 6.

For the right part of Figure 4 we have:

sin θ r 2 = AC R (1) sin θ r 2 = GD h (2) cos θ r 2 = GF h (3)
This is used to compute the angle of the incident ray:

tan (i) = AC EC = R sin θ r 2 h -R 1 -cos θ r 2 (4)
Meanwhile, we have:

ψ = π 2 - ε -θ r 2 ⇒ sin (ψ) = cos ε -θ r 2 (5)
Therefore, using Snell's law, the velocity ratio can be expressed as:

V e V d = cos ε -θ r 2 1 sin arctan R sin ( r 2R ) h-R(1-cos( r 2R )) ⇒ arcsin V d V e cos ε -θ r 2 = arctan R sin r 2R h -R 1 -cos r 2R ⇒ R sin r 2R h -R cos r 2R = tan arcsin V d V e cos ε 2 - r 2R
which gives:

h = R sin r 2R tan arcsin V d V e cos ε 2 -r 2R + R 1 -cos r 2R (6)
A consequence of this equation is that for a given epicentral distance, there exists a unique pair of range and diffusive layer thickness for a given velocity ratio (i.e., a ray parameter that satisfies both). In Figure 5 we show the values obtained for the range of the seismic ray in a unit layer thickness (h = 1) and epicentral distance ε = 20 • , 50 • and 100 • as a function of the velocity ratio,

V d /V e .

Computation of energy envelopes

A common approach in seismic coda analysis is to model the energy envelope of the seismic wavefield, which discards information about polarity and phase in favor of the shape of the coda decay. To fit the shape of the S-coda energy envelope, we use the [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF] equation to compute a theoretical energy envelope of the S-coda waves, considering an impulse at the source:

i (t) = 4 πξ H th exp - r 2 ξ H t - wt Q ∞ ∑ n=1 α n 2α n + sin 2α n exp - tξ V α 2 n 4h 2 (7)
where ξ H and ξ V are the horizontal and vertical diffusivity components, t the time, h the diffusive layer thickness, r the seismic ray range in cylindrical coordinates, w the frequency calculated as w = ( f min + f max )

2

, where f min and f max are respectively the minimum and maximum frequency of the selected envelope from the data, Q the attenuation factor, and α the positive roots of the equation:

α tan α = 4hυ ξ V ( 8
)
where υ is the seismic velocity in the underlying elastic halfspace. This means that υ corresponds to the V e term of Equation 6, so that holding fixed the velocity in the diffusive layer, υ will increase with the diffusive layer thickness, h. Therefore, the boundary condition that is described by Equation 8 depends on the diffusive layer thickness.

This equation was developed for lunar impacts and the shown form is valid only for shallow seismic events, with depth z = 0 in cylindrical coordinates. As we discuss later, this assumption makes the methodology less appropriate when applied analyzing coda of deep marsquakes.

Diffusivity computation and its dependence on the geometry

We assume a diffusivity in the diffusive layer D = ξ H = ξ V although we note that in practice the ratio of the horizontal and vertical diffusivity is typically greater than one owing to the additional seismic energy contributed by surface waves to the horizontal component of motion. Dainty et al. (1974a) defined the relationship of the diffusivity, D, with the free mean path, l, and the velocity of wave propagation, υ, as following:

D = υl 3 (9)
The transport mean path, l is given by:

l = (σ n) -1 ( 10
)
where σ is the cross-sectional area and n the number of particles per unit volume. If we consider two layers with the same seismic wave velocity, υ 1 = υ 2 , we find that the diffusivities D 1 and D 2 are proportional to l. In order to compute its value we can consider a cuboid of dimensions h (the diffusive layer thickness) and r (the range of the seismic ray). The number of events corresponding to the cross section of this cuboid is equal in all directions, as we assumed earlier that ξ H = ξ V . Given that n is the number of particles per unit volume, corresponding to a cross sectional area σ , the number of scattering events, N, in a cross-sectional area S = hr is:

N = Sn = hrn (11)
Starting with an impulsive signal, in order to obtain the same envelope for υ 1 = υ 2 but different size of the diffusive layer, we need to have N 1 = N 2 . Using the Equation 11, we have:

N 1 = N 2 ⇒ h 1 r 1 n 1 = h 2 r 2 n 2 ⇒ n 1 n 2 = h 2 r 2 h 1 r 1 (12)
Then we solve this on the basis of the definition of diffusivity (Equation 9) to obtain:

D 1 D 2 = υl 1 υl 2 = l 1 l 2 = σ n 2 σ n 1 = n 2 n 1 (13) D 1 D 2 = h 1 r 1 h 2 r 2 ⇒ D 1 = h 1 r 1 h 2 r 2 D 2 (14)
This relationship between the thickness of the diffusive layer and therefore the range (which is itself uniquely related to the velocity ratio) of the seismic ray, leads to a tradeoff with the diffusivity that is discussed in the results. Thus, in order to explore the effect of each parameter, we perform a grid search over the model space. Similarly, Lognonné et al. ( 2020) used 3 different crustal thicknesses (h = 20, 40 and 60 km) to investigate the respective diffusivity for the Martian crust.

In Figure 6 we show the effect of this tradeoff on the computed energy envelopes. On both left and right panels we use the same velocity ratio between the diffusive and elastic layer, the same Q attenuation factor and the same frequency, w, as these variables are contained in Equation 7. We compute energy envelopes for 3 different layer thicknesses, h = 10, 20 and 30 km and solve Equation 6 for the self-consistent value of r. In the left panel, we fix the diffusivity to d = 0.1 km 2 /s, and can observe that the coda decay duration is longer for a thicker diffusive layer.

In the right panel, we adjust the diffusivity for each layer, using Equation 14 and the thickness of 20 km as a reference. We observe that for all three layer thicknesses, the adjusted diffusivities produce identical envelope shapes. Due to this complete tradeoff between diffusivity and thickness, we only need to compute the results for a given layer thickness, obtain the other parameters of Equation 7 for that layer thickness, and then adjust their values accordingly for the thickness of the diffusive layer.

Reverberation dependence

A possible complication to our average ray path assumption is the effect of rays that do not follow the direct path of the refracted ray, but rather reverberate within the scattering layer(s). As described earlier in this section, with the approximation of the average seismic ray, due to the reflection on the interface with the elastic layer, the reverberations in the diffusive layer can be considered part of the diffusion of the seismic waves. Therefore, a hypothetical n number of reverberations of range r can be modeled as a unique seismic ray of range n • r in one layer of a given thickness. The diffusivity in this case should be adjusted with the use of Equation 14for the new range of the seismic ray.

In order to observe the effect of this adjustment, we perform a test that is shown in Figure 7.

On the left panel, we show in blue the envelope for a single ascending ray in a 20 km layer, with a Q = 1000 and d = 1 km 2 /s. Using the same layer thickness and Q, we compute the energy envelope for a seismic ray that is generated on top of the diffusive layer and is reflected on the interface with the elastic layer. If we define the range of the single ascending ray as r, in the case of a single reflection, we are computing two distinctive envelopes. The first is the result of an impulse, using the Equation 7, a range r/2 and diffusivity 2 km 2 /s (with the use of Equation 14).

The second is the result of the same computation, however instead of inputting an impulse, we use the result of the first envelope computation. The result of this single reflected ray is shown in red.

We apply the same methodology accordingly to the cases of a double and triple reflection in the diffusive layer, showing the respective results in yellow and purple. On the right panel, the same computations are performed for an initial diffusivity d = 5 km 2 /s.

We observe that the computed envelopes do not differ significantly by adding more reverbera-tions to the ray path. More precisely, the computed deviation of every case of reverberations versus the case of the ascending ray is shown in Table 2. This deviation does not appear to depend on the number of the reverberations or the diffusivity, whereas in all cases it is much less than 0.1 which is approximately the best misfit that we find in the results presented in the next section.

It is important to note that there is no ballistic wave in the diffusion model, hence no reflected waves in the traditional sense. As illustrated by [START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda q[END_REF], the range of validity of the diffusion models extends to a point where the mean free path, l, is longer than the thickness of the scattering/diffusive layer; modeling beyond this range of validity is outside the scope of this study.

RESULTS

Here we examine the fit between the observed spectral energy envelopes from the Martian data and the theoretical scattering model. We perform a grid search by using a range of values for the parameters of Equation 7, which define the characteristic diffusivity and scattering attenuation in the shallow Martian lithosphere. Models that minimize misfit with the data are considered possible structures for the interior of Mars.

Previous studies [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF] have noted that it is impossible to constrain the diffusive layer thickness based on the S-coda wave analysis alone, because layer thickness has a tradeoff with the diffusivity. As discussed in the Methodology section, we do not need to vary the layer thickness in the computations because the results for any desired layer thickness can be calculated given the range of the seismic ray and the diffusivity. In our grid search, we use a fixed layer thickness, h = 20 km, and adjust the other parameters (thickness, range) accordingly. Figure 8 shows the necessary adjustment of the diffusivity, depending on the layer thickness.

The range of the seismic ray depends on the epicentral distance and the velocity ratio between the diffusive and elastic layer, and is obtained by solving Equation 6 for r. The epicentral distance for all the studied events is given in the Seismic Catalog (InSight Marsquake Service, 2021) and reproduced in Table 1, and we perform our investigation for 5 different velocity ratios, V d /V e = 0.15, 0.18, 0.20, 0.25, 0.30. One can note that the value for this ratio is much smaller than the velocity ratio expected between typical seismic discontinuities; for example, the crust and the mantle on Earth has a velocity ratio of around 0.6, and other crustal layers may be even higher.

However, because we base our modeling approach on the average seismic ray through the nonscattering and scattering medium, we expect the apparent velocity in the diffusive layer to be lower due to the scattering. Finally, we search in the range of Q s = 100 to 2000 for the quality factor of in the scattering layer.

For each case of the aforementioned parameters, we select a time window for the computation of the misfit between the data (spectral envelopes) and the prediction (computed energy envelopes).

To calculate misfit, we first align the observed and modeled envelopes on their peak, and trim the model data series in the appropriate time window that corresponds to the S-wave arrival and coda of the data, as explained in the Data Processing section. We finally compute the Root Mean Square Error (RMSE) between the model and data:

RMSE = ∑ N i=1 (d i -s i ) 2 N ( 15
)
where d i is the observed spectral envelope and s i the synthetic one.

In Figure 9 we show the data and computed models for an event of each frequency type. The black curve shows the spectral envelope computed from the vertical component of the velocity seismogram. The red area indicates the margins of normalized amplitudes for a range of models that provide an RMSE lower than the indicated threshold on the top of each example. It is noted that the RMSEs are much higher for the LF and BB events, whereas the HF and VF events show a better match to the model prediction. This is consistent with the suggested location and focal depth of the events, as was analyzed by [START_REF] Giardini | The seismicity of mars[END_REF]. The LF and BB events are located at teleseismic distances and are inferred to be deep, possibly sub-crustal marsquakes. Their ray paths would therefore travel longer in the elastic region of the Martian interior than the respective HF and VF events. The latter, which are typically located at shorter distances and assumed to be events near the surface, would have waves that are propagated through the highly diffusive layer (or region) near the Martian surface.

The grid search results for every event are examined in a summary plot, and an example is shown in Figure 10 for the event S0231b (HF). The results are shown for a layer thickness of h = 20 km. The diffusivities shown in this summary should be adjusted for each desired layer thickness, as described earlier in this section and shown in Figure 8. The white curves show the lowest misfit for every pair of Q and diffusivity values and correspond to the curves that are shown on the right bottom side, where the RMSE is plotted as a function of Q for each case of velocity ratio. The gray dashed line in this subplot corresponds to the best misfit among all the velocity ratios tested.

We observe that for higher velocity ratio, there is a narrow region of low misfit, which indicates a preference for a specific small range of Q. For smaller values of the velocity ratio, this range of preferred Q increases and the associated curve reaches a flat region for the lower RMSEs. The shape of the gray dashed line shows that we cannot choose a specific attenuation factor based on this analysis. Importantly, its shape is not identical for every event family (LF, BB, HF and VF)

which allows us to deduce information about the properties of Mars by a comparative analysis.

The results of that analysis are shown in the bottom right part of Figure 10 for each event, and are presented collectively, for the events of each type in Figure 11. The best and worst misfit (lower and higher RMSEs) as a function of Q are shown in red, whereas the mean value of all the curves is shown in blue.

Due to the small number of lower-frequency type events (LF and BB), the results for these event types can be considered as more uncertain. For the LF events, the lowest misfit corresponds to the curve for event S0409d and the highest misfit for event S0189a. As seen in Table 1, we do not observe any correlation of the envelope fit and the frequency content of the data envelopes or the epicentral distance of the events. For the BB events, there are only 2 computed curves, with the minimum misfit corresponding to event S0235b and the maximum to event S0185a. This could be evidence that the misfit, in the flat region of higher Q, correlates with the epicentral distance of the events, but the same correlation is not found across all the event types, and remains unconstrained. For the LF and BB events, we note the inability of our modelling approach to provide good envelope fits independently of the parameter ranges used in the grid search.

On the other hand, more information about the structural properties of Mars is provided through the analysis of the HF and VF events. As shown in Figure 10 for the HF event S0231b, at lower values of Q, the spectral envelopes can be fit with a smaller diffusivity, while at higher Q values, the diffusivity must also increase; best fits (lowest RMSE) are found for intermediate Q values.

This local minimum in the RMSE as a function of Q is even more apparent for VF events, where optimal fits are provided by Q values in the 400 -640 range. For VF events, we also find that as the frequency content of the events increases, a greater quality factor provides better fits to the data. However, this is based on the analysis of only 5 events and more event data are necessary in order to draw any firm conclusions from the frequency and scattering quality factor correlation.

The same correlation is not observed in analysis of the HF events.

Due to the tradeoff between the layer thickness and the diffusivity, which is discussed in the Methodology section and summarized by Equation 14, it is not possible to constrain a specific diffusivity in the elastic layer. This tradeoff is further illustrated in Figure 12 using the results of the analysis of VF event S0128a. On the left panel, the minimum RMSE is plotted as a function of Q. The different color curves correspond to different velocity ratios, which result in a different range for the seismic ray and therefore distinct ray lengths in the diffusive layer, as indicated in the legend. We observe no preference for any specific diffusive layer thickness as they can all satisfy the envelope equally well, for a different choice of diffusivity. These corresponding diffusivities are shown in the center panel, with increasing values for increasing Q. This feature is a direct consequence of Equation 7. To fit the data envelope, there is a large range of acceptable diffusivity values that trade off directly with the chosen layer thickness, as shown on the right panel of Figure 12.

The complete trade off between layer thickness and diffusivity is a major obstacle for the interpretation of our results in terms of scattering layer thickness and strength. The obtained diffusivity results that correspond to a specific diffusive layer thickness can be adjusted at will, by following the relationship of Equation 14, as it is demonstrated by the test shown in Figure 6. However there is another element of the analysis that can be used to constrain the structure of the diffusive part of the Martian lithosphere. As shown in Figure 12, the diffusivity varies with the change of the velocity ratio between the studied diffusive layer and the underlying elastic half-space. When V d /V e increases, the diffusivity will increase as well. In order to investigate if this is another artifact due to the tradeoff between the diffusivity and the dimensions of the seismic ray path in the diffusive layer (defined by h and r) we use the results of the VF events, as shown in Figure 12 for event S0128,a and test if the computed diffusivities can be obtained by only using the Equation 14 for a constant layer thickness and the respective range of the seismic ray, corresponding to different velocity ratios. The computed diffusivities depend on the change of the range, r, of seismic ray, however not linearly but quadratically, which means that for any given diffusivity D 1 , for a velocity ratio (V d /V e ) 1 with corresponding range r 1 and a layer thickness h, there is a diffusivity D 2 for (V d /V e ) 2 that gives a seismic ray range r 2 for the same layer thickness and their relationship is:

D 1 D 2 = r 2 1 r 2 2 (16)
This relationship can be obtained through the joint solution of Equations 6 and 14. Therefore, for the interpretation of the data, we need to analyze the range of best fitting diffusivities and Q pairs for one given velocity ratio and layer thickness and therefore range of the seismic ray.

DISCUSSION

In our analysis we investigate the scattering properties of the Martian interior, based on the computation of energy envelopes used previously for lunar impacts by [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF]. By systematically investigating the effects of all the model parameters (7), we establish the existence of key tradeoffs between these parameters. More precisely, we show that there are trade offs between the dependence of the diffusivity, the diffusive layer thickness, the velocity ratio between the diffusive layer and the underlying elastic half space, and the range of the seismic ray in the diffusive layer.

This means that by knowing any one of these parameters independently, we can use our modeling to place definite bounds on parameters controlling the scattering in the Martian interior. On the other hand, if we do not have independent constraints for these parameters, we end up with a multidimensional space of possibilities for models that fit the data. [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF] showed results for specific models of scattering, with separate analyses assuming a diffusive layer thickness of h = 20, 40 and 60 km. More precisely, they performed a preliminary analysis of VF event S0128a, LF event S0173a and BB event S0235b. The key purpose of their work was to demonstrate the compatibility of observed envelope shapes with a multiple-scattering origin. In the case of the VF event, they used forward modeling based on radiative transfer equations in a few sets of statistically uniform random models with ad-hoc statistical properties. These authors inferred that a diffusivity of the order of D = 90 km 2 /s at a frequency of 7.5 Hz was compatible with the observations. However since the modeling did not include any depth dependence, our work shows that this value may be an overestimate. In the case of the LF/BB events, Lognonné et al. ( 2020) considered a simplified model where a single plane wave impinges vertically on a scattering crust from below. While their approach bears some similarity with the present one, there are some important differences with our work: the bottom of the scattering layer in [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF] coincides exactly with the Moho, whereas any scattering effect on the downgoing part of the ray is neglected. In that study it is estimated a broad range of values for the diffusivity at 0.5 Hz, from 200 km 2 /s to 2000 km 2 /s. In the present work, we considerably expand the initial dataset, thereby covering a broader range of frequencies and epicentral distances. We also extend the range of values of the parameters that control scattering. Our study investigates in details the possible trade-offs between assumptions for scattering strength, crustal thickness, velocity contrast at the Moho and absorption. Thereby, we offer a more comprehensive view of the current uncertainties on the scattering properties in the Martian lithosphere.

Through our analysis of envelope shapes, we find that many models could fit the data, and that drawing conclusions based on a single unique model that can fit the data would be misleading, as it would be reflect assumption(s) made to remove tradeoffs among several parameters that control scattering. For example, previous numerical modeling of scattering proposed candidate structures that fit the InSight data [START_REF] Van Driel | High frequency seismic events on mars observed by insight[END_REF], for a specific layer thickness of h = 10 km, velocity ratios v d /v e > 0.5 and scattering ranges from 10% to 100%. These ranges of parameters are all found to yield acceptable fits to the coda envelopes analyzed in our study. More precisely, if we apply Equations 6 and 16 to the results that we present in Figure 12 for HF event S0128a, we find that for this range of parameters, the diffusivity that provides the best fit is D = 0.7 km 2 /s, which is in agreement with the findings of the van Driel et al. ( 2021) paper.

In Figure 13, we present a comparison between envelopes computed through 2D numerical wave propagation simulations of van Driel et al. ( 2021) and those computed in our study using the analytical methodology based on the theory by [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF]. For the estimated range of diffusivity used in that study (D = 0.5 -0.7 km 2 /s), our predicted coda decays provide an imperfect but good fit to their synthetic coda envelope, as seen in the left panel of Figure 13. In the right panel, we increase the range of diffusivity values to D = 0.3 -0.7 km 2 /s, and plot the coda decays predicted by the approach developed in our paper. We find that the lower diffusivity values provide improved fits to the S-coda envelopes computed through numerical wave propagation simulations, mainly for epicentral distances greater than 15 • (receiver index greater than 7 in Figure 13). At smaller epicentral distances, higher diffusivity values fit the modeled data better, which may be the reason that we find a good agreement for the results for event S0128a, a Very High Frequency event at an estimated epicentral distance of 7.79 • (in the middle range between receiver index 3 and 7). This comparison with coda decays obtained through numerical wave propagation simulations of [START_REF] Van Driel | High frequency seismic events on mars observed by insight[END_REF] demonstrates that the relationships developed in this study yield correct estimates for the parameters that control scattering.

In their study of the lunar interior, [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF], who originally used the modeling equation that was chosen for our analysis, suggest a unique model that can fit the ensemble of the data. However, in their model they suggest an apparent thickness for the diffusive layer which varies with the frequency of the examined data. In order to do this, they use a factor for the amplitude of the body waves generated by the Lunar Module and the Saturn-IV B impacts, that was computed by the study of [START_REF] Toksöz | Velocity structure and properties of the Lunar Crust[END_REF]. We showed that modeling based on Equation 7yielded better fits to the observed S-coda evelopes of higher frequency Martian events (HF and VF). This finding is consistent with previous studies that suggested that these are shallow events [START_REF] Giardini | The seismicity of mars[END_REF][START_REF] Van Driel | High frequency seismic events on mars observed by insight[END_REF]. We also showed that the approach based on the mean ray path can describe with fairly small errors eventual reverberations of the diffusive waves in a shallow diffusive layer. As shown in Table 2 and Figure 7, these errors are much smaller than the RMSE between the computed energy envelopes and the spectral envelopes of the data.

In this study we examine only the first arriving S-wave signals in our data but excluded the Pwave coda which also contains complementary information on scattering. The P-waves propagate with a different velocity and frequency content, and future analysis of their coda properties could provide an independent constraint on scattering in Mars. Furthermore, we normalize the maximum amplitudes of each event envelope to unity, removing information on the absolute energy loss.

While the magnitude of the energy loss in the elastic layer does not strongly affect the envelope shape, an analysis of absolute amplitude at varying event distances would provide constraints on the intrinsic Q of the mantle, a task we do not explore further here. Finally, we assumed the elastic region underlying the scattering layer to be a homogeneous half-space, which does not affect the envelope shapes and therefore the coda decay. This means that we do not compute the actual seismic rays but only an average interpretation of their paths. Future work could combine data analysis to constrain the model space of scattering properties, with more sophisticated full wave propagation modeling to further refine the fits to marsquake waveforms.

We can interpret the success of modeling VF and HF events and the relative inability to model LF and BB events in terms of the likely source characteristics for these event types [START_REF] Giardini | The seismicity of mars[END_REF]. Compared to the HF or VF envelopes, we observe a very rapid coda decay for the LF events and slightly longer, but still rapid decay for the BB events. Our approach for the computation of the envelopes in a diffusive layer cannot be effective when we try to model the body waves that propagate in the elastic part, which would be expected for body waves from deep marsquakes.

For deep events, multiple scattering that can be modeled as a diffusion process happens only in the vicinity of the station, when the impulse of seismic energy broadened only by the effects of attenuation arrives from below. Thus, we interpret the inability of our model to fit the coda of LF and BB events to be further evidence that these events are deep marsquakes.

This might be due to the fact that the waves travel shorter distances in the diffusive layer when they occur deeper in the Martian interior, which has also a correlation with their frequency content. An argument that supports this hypothesis is the better fit of the equation developed for lunar impacts [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF]) (i.e. shallow and surficial) to higher frequency events. It is in coherence with the suggestion that lower frequency events occur in greater depths. However, the modeling of these events is not able to show a critical result to constrain the thickness of the diffusive layer. In addition, regarding the poor fit of LF and BB events, we know now that they are composed of multiple energy injections that arrive in the coda. These injections are not taken into account in the modeling and this contributes to the difficulty fitting the data. It is additionally worthy to note that [START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF] performed an analysis of the LF and BB events, showing that respective spectral envelopes can be modeled through the analysis of the ratio between the ballistic and S-coda waves, as well as the coda decay, but they do not follow the diffusion model.

The epicentral distances of the HF events have a small range in variation, between 21.4 • and 28.4 • (see Table 1) which limits our ability to study how scattering changes over distance. Translated into km, these values correspond to a distance between 1262.6 km and 1675.6 km. This means that the recorded waves should either cross the crust-mantle boundary and travel in an elastic regime in the lithosphere, or become trapped and reverberate in the crust, as suggested by [START_REF] Giardini | The seismicity of mars[END_REF] and [START_REF] Van Driel | High frequency seismic events on mars observed by insight[END_REF]. In the case of many reverberations, the apparent speed in the diffusive layer should be very high according to our analysis. This high seismic velocity corresponds to a very low Q in our results and very low diffusivity, which corresponds to a region where our results are saturated. However, for the VF event S0128a, located at a relatively small epicentral distance (7.79 • ), we find that a low Q = 200 -300 fits the data better than a higher Q. It is therefore unclear if the crustal waveguide, or properties of the scattering in the crust are producing these difference between more distant or closer marsquakes. As more HF events are recorded at different distances, this behavior can be investigated in more detail.

The analysis of the VF events shows that they are better modeled with Equation 7, as would be expected if their sources were indeed shallow or surficial [START_REF] Giardini | The seismicity of mars[END_REF]. Epicentral distances of VF events vary widely from 6.44 • to 36.8 • . Despite this distribution in terms of distance, we do not observe a correlation between their distance and their S-coda decay time or for inferred values of diffusivity, as would be expected for a near surface layer of some given thickness and scattering properties. One possibility is that the VF events have a distribution of azimuths with respect to the InSight lander position and thereby sample very different scattering structures.

Unfortunately, for most of the VF events, it is not possible to robustly determine backazimuth.

Nevertheless, future determinations of VF backazimuths would make it possible to infer lateral variations in the Martian crust, and argue against a single, uniform diffusive layer.

To interpret the range of structures that fit the obtained results, we can assume for example, a relatively thin diffusive layer confined to the shallow-most crust. Underneath, the elastic halfspace contains part of the diffusive region of the Martian crust and upper mantle's structure. The diffusivity in the thinner top layer should be smaller than the one obtained for a thicker one. This tradeoff suggests that the number of scattering events depends not only on the number density of scatterers in a given cross section (Equation 14) but also depends on the square of the range of the seismic ray (Equation 16).

Given this interpretation, a last question is whether the examined diffusive layer structure, with its lateral variations, corresponds to a part of the Martian crust or extends deeper in the crustmantle or even in the upper mantle. If we assume that future analyses or events yield reliable backazimuth estimates, the answer to this question depends on the level of the general knowledge for Mars interior and more precisely the structure of the lithosphere, as this approach will be able only to constrain the thickness of diffusive layer but not the thickness of the crust, which is a matter of debate in the literature. [START_REF] Wieczorek | Global character of the martian crust as revealed by insight seismic data[END_REF] performed a review of the studies that defined the average crustal thickness of Mars. The suggested values vary from H C < 29 km for a model of isostatically compensated crust in Hellas Planitia [START_REF] Wieczorek | Thickness of the martian crust: Improved constraints from geoid-to-topography ratios[END_REF], to H C < 115 km for a model with viscous relaxation of dichotomy boundary and Hellas basin [START_REF] Nimmo | Estimates of martian crustal thickness from viscous relaxation of topography[END_REF]. Moreover, they suggested a crustal thickness in the vicinity of InSight either 20 km or 37 km through the assumption of a 2-layered or 3-layered crust model respectively.

More recently, Knapmeyer-Endrun et al. ( 2021) computed autocorrelations and receiver functions using InSight data and suggested an average crust thickness varying between 24 -70 km, with the presence of either one or two seismic interfaces in the Martian crust. The relevance of these crustal thickness constraints for the interpretation of results presented in our study hinges on whether the seismic waves of the HF and VF events are crossing the crust-mantle interface or if they travel only within the diffusive layer considered as part of the crust defined by intra-crustal interfaces.

Information about the energy loss towards the inner depths of the planet, which can be provided through analogies of the expected amplitudes as it was done in previous works (for example [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF] used a known analogy for the amplitudes of the Lunar Module and Saturn IV artificial impacts, provided by [START_REF] Toksöz | Velocity structure and properties of the Lunar Crust[END_REF] in order to suggest a unique structure model) can be a valuable element to constrain this feature. Furthermore, it will be useful if future experiments are performed in a region close to the InSight seismic experiment (more precisely in a range of around 30 • , as this is an average distance of the High Frequency events) and this geographical setting will allow a joint analysis and further interpretation of the currently available InSight seismic data. The existence of such a network will improve the ability of phase peaking and location identification of the events and therefore it will give an extra constrain for an analysis similar to this study's, which is the structure of each event's waves propagation, with more data coming from events that are now characterized of lower quality in the Seismic Catalog [START_REF] Clinton | The marsquake catalogue from insight, sols 0-478[END_REF].

Mars appears to be intermediate between the Earth and Moon in terms of seismic scattering and attenuation, since Martian seismograms exhibit a shorter coda durations than lunar seismograms [START_REF] Van Driel | High frequency seismic events on mars observed by insight[END_REF]. However, the HF and VF Martian events have long codas and exhibit some resemblance to Moonquakes and strongly scattered Earth seismograms, while the LF and BB events resemble regional tectonic events on Earth. Like the Moon, the origin of the scattering on Mars likely lies in the crust or uppermost lithosphere. On the Moon, the scattering is produced by impact processes that have produced a shallow layer of regolith and deeper megaregolith of highly fractured bedrock. The scattering properties of lunar events were measured by [START_REF] Gillet | Scattering attenuation profile of the Moon: Implications for shallow moonquakes and the structure of the megaregolith[END_REF], who derived a model of the scattering and attenuation properties of the Moon using diffusion theory. They found very low wave diffusivity (D ≈ 2 km 2 /s) in the uppermost 10 km of the Moon. They noted that these values correspond to some volcanic areas on Earth, which are the most heterogeneous regions on our planet. Below the surface layer, the diffusivity rises slowly up to a depth of 80 km, where it increases abruptly by about one order of magnitude. [START_REF] Gillet | Scattering attenuation profile of the Moon: Implications for shallow moonquakes and the structure of the megaregolith[END_REF] suggested that the megaregolith corresponds to the region of low diffusivity, and that it is 100 km thick (much larger than previous estimates). When looking at the seismic layers at Mars, there is a low velocity surface layer [START_REF] Knapmeyer-Endrun | Thickness and structure of the martian crust from insight seismic data[END_REF][START_REF] Lognonné | Constraints on the shallow elastic and anelastic structure of mars from insight seismic data[END_REF], which probably represents the ejecta rubble and severely cracked rock produced by lunar meteorite bombardment [START_REF] Goins | Lunar Seismology: The Internal Structure of the Moon[END_REF]. Our Martian data require low diffusivity (D is generally lower than 1.5 km 2 /s), which suggests that the diffusivity in the top 10 km beneath Elysium Planitia on Mars are similar to the low diffusivity found on the Moon. However, there is weak geological evidence on Mars for a thick megaregolith layer, which is further substantiated by the existence of the Low Frequency and Broadband Marsquake events [START_REF] Giardini | The seismicity of mars[END_REF] that appear to occur below the scattering layer. Therefore, Mars appears to be more complicated than a simple intermediate between Earth and the Moon. Instead, it shares some of the properties of these two bodies, the only other seismically-investigated bodies in the Solar System.

CONCLUSION

We investigated the seismic attenuation in the Martian crust and upper mantle by examining the Swave codas of a series of InSight detected marsquakes. For our investigation, we used the spectral envelopes 21 marsquakes in 4 different event families -classified by their frequency content -with source parameters from the Seismic Catalog (InSight Marsquake Service, 2021). We assumed a diffusive layer over an elastic half-space model, and computed the mean raypath of the seismic waves from the shallow source to the station for a given epicentral distance and for free variables of scattering, Q, diffusivity, and velocity ratio between the diffusive and the elastic layer.

In our study, we observed that the Low Frequency (LF) and Broadband (BB) events, with frequency content below the threshold of the 1 Hz tick noise, could not be fit by our model. The spectral envelopes of the S wave codas of these events showed a very rapid decay which suggests that they do not have an extensive propagation path in the diffusive layer. This observation is in agreement with the suggestion of previous studies [START_REF] Giardini | The seismicity of mars[END_REF]) that these events are deep marsquakes and travel through the upper mantle of Mars.

Based on the results of the High Frequency (HF) and Very High Frequency (VF) events, we observed a range of possible paths and diffusivities that can satisfy the data, and we investigated the tradeoffs between the parameters in the modeling equation [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF]) that controls the shape of the energy envelope for the events. The analysis of these tradeoffs shows that the S-wave coda do not uniquely constrain the depth of the diffusive region in the Martian crust and the upper mantle. Our analysis of HF and VF events, is consistent with previous studies [START_REF] Giardini | The seismicity of mars[END_REF][START_REF] Van Driel | High frequency seismic events on mars observed by insight[END_REF] which argued that that these events are shallow-sourced.

However, the observation that the lower frequency event families cannot satisfy the model, showing a very rapid S-coda decay, suggests the possibility that the Martian lithosphere may differ compared to the lunar or terrestrial one. The diffusive region on Mars is comparable to the lunar regolith, however the regolith on Mars is not extended to great depths, as demonstrated by deep marsquakes that appear to propagate in the elastic region of the Martian lithosphere. We find that the equation used by [START_REF] Dainty | Seismic scattering and shallow structure of the moon in oceanus procellarum[END_REF] to fit moonquakes is not able to fit the Martian data at all the frequency ranges. Therefore, we deduce that the scattering structure of the Martian lithosphere is not similar to the Moon, and that Mars is unlikely to have a deep megaregolith.

The results of this study illustrate the challenges of working with single station seismic data where independently determined event location information, including distance, azimuth, and depth are crucial for understanding the lateral variation in seismic properties of a planet.

DATA AND RESOURCES

Seismic data used for this study were collected as part of the Seismic Experiment of Internal Structure (SEIS) [START_REF] Lognonné | Seis: Insight's seismic experiment for internal structure of mars[END_REF] of the NASA InSight Mission to Mars [START_REF] Banerdt | Insight: a discovery mission to explore the interior of mars[END_REF].

They can be obtained from the IRIS Data Management Center (https://www.iris.edu/hq/sis/insight), the NASA PDS Geoscience Node (InSight SEIS Science Team, 2019) and the IPGP SEIS Data portal (InSight Mars SEIS Data Service., 2019).

Ludovic Margerin

Institut The mean ray path (red) of the seismic waves from the source (A) to the station (D) through an elastic layer that overlies a diffusive layer. The refracted waves case is shown on the left and the reflected waves on the right. The thickness of the diffusive layer is noted with h, the range of the seismic ray in the diffusive layer with r in length units and θ r in radians, the S-wave velocity in the diffusive layer with V d and in the elastic layer with V e . The epicentral distance is noted with ∆ in length units and ε in radians. R is the planet's radius.

Figure 5

The range of the seismic ray in the diffusive layer, for a layer thickness h = 1 and epicentral distances of 20 • , 50 • and 100 • for a range of velocity ratio between the diffusive and elastic layer from 0.1 to 1.

Figure 6

The energy envelopes obtained for an impulse source at the surface. In both panels, the blue, red and yellow curves correspond to layer thicknesses of 10, 20 and 30 km respectively. On the left panel, the diffusivity is D = 0.1 km 2 /s for all the examples. We observe longer decay rates for bigger diffusive layer thickness. On the right panel, the diffusivity is adjusted using Equation 14. With this adjustment, which follows the tradeoff between layer thickness and diffusivity, the computed envelopes remain unchanged.

Figure 7

The effects of adding multiple reverberations into the envelope calculation. The blue line represents a seismic ray traveling from the discontinuity between the diffusive and the elastic layer, upwards towards the station. The red line represents a seismic ray of a surface event with one reverberation and the yellow and purple line represent 2 and 3 reverberations respectively. In both panels, the diffusive layer thickness is set to h = 20 km and attenuation factor Q = 1000. The diffusivity on the left panel is D = 1 km 2 /s and in the right panel D = 5 km 2 /s.

Figure 8

The tradeoff between diffusivity and layer thickness. Using a diffusivity D = 1 km 2 /s for a thickness of h = 20 km, we use Equation 14 in order to show how the diffusivity should be adjusted, depending on the layer thickness, in order to obtain the same results in our modeling.

Figure 9

The fit of the scattering diffusive model to spectral envelopes of the data. The spectral envelopes of 4 seismic events, one for each family, are shown in black. The synthetic envelopes that provide an RMSE lower than a specific threshold (0.4 for the LF, 0.16 for the BB, 0.14 for the HF and 0.1 for the VF events) are shown in red. We observe that the modeling approach works better for the HF and VF events, which are considered to be located at closer epicentral distances and be sourced near the surface.

Figure 10

Results for models fits of HF event S0231b. The 5 colormaps show the Root Mean Square between the data spectral envelope and the computed modeled envelope for the velocity ratio shown on the top of the subplot, the Q on the y-axis and the diffusivity given in km 2 /s for a layer thickness of h = 20 km. The white curves note the best misfit for every pair of Q and diffusivity. The best misfit is also shown on the right bottom side, in respect to Q for each case of velocity ratio. The gray dashed line in this subplot corresponds to the best misfit among all the velocity ratio associated curves.

Figure 11

The curves of the best misfit as a function of Q for each event family. The dashed lines correspond to the results of every event, whereas the minimum and maximum values of the former are shown in red and their average value is shown in blue. The specified frequency range corresponds to the lower and upper frequency that was used for the filtering of the ensemble of the data of each event family.

Figure 12

The analysis of the results for VF event S0128a. a) The best fit between the data and the computed envelopes for each Q is shown. Different colors correspond to the velocity ratio between the diffusive and elastic layers, which controls the length of the ray path in the diffusive layer, as indicated in the legend. The green cross shows the minimum RMSE for a Q = 400. b) The corresponding diffusivities are shown and we observe an increasing diffusivity that satisfies the Figure 6: The energy envelopes obtained for an impulse source at the surface. In both panels, the blue, red and yellow curves correspond to layer thicknesses of 10, 20 and 30 km respectively. On the left panel, the diffusivity is D = 0.1 km 2 /s for all the examples. We observe longer decay rates for bigger diffusive layer thickness. On the right panel, the diffusivity is adjusted using Equation 14. With this adjustment, which follows the tradeoff between layer thickness and diffusivity, the computed envelopes remain unchanged. that provide an RMSE lower than a specific threshold (0.4 for the LF, 0.16 for the BB, 0.14 for the HF and 0.1 for the VF events) are shown in red. We observe that the modeling approach works better for the HF and VF events, which are considered to be located at closer epicentral distances and be sourced near the surface.

  Figure 4

Figure 1 :

 1 Figure 1: The spectral envelopes of example events for each of four different event types defined

Figure 2 :

 2 Figure 2: The complete dataset of seismic events used in this study for vertical velocity seis-

Figure 3 :

 3 Figure 3: Schematic showing the approximate average seismic raypath from the source to the

Figure 4 :

 4 Figure 4: The mean ray path (red) of the seismic waves from the source (A) to the station (D)

Figure 5 :

 5 Figure 5: The range of the seismic ray in the diffusive layer, for a layer thickness h = 1 and

Figure 7 :

 7 Figure 7: The effects of adding multiple reverberations into the envelope calculation. The blue

Figure 8 :

 8 Figure 8: The tradeoff between diffusivity and layer thickness. Using a diffusivity D = 1 km 2 /s

Figure 9 :

 9 Figure 9: The fit of the scattering diffusive model to spectral envelopes of the data. The spectral

Figure 10 :

 10 Figure 10: Results for models fits of HF event S0231b. The 5 colormaps show the Root Mean

Figure 11 :

 11 Figure 11: The curves of the best misfit as a function of Q for each event family. The dashed

Figure 13 :

 13 Figure 13: The envelopes obtained with the analytical method in this study (red, orange) are

ACKNOWLEDGEMENTS

We acknowledge NASA, CNES, their partner agencies and Institutions (UKSA, SSO, DLR, JPL, IPGP-CNRS, ETHZ, IC, MPS-MPG) and the flight operations team at JPL, SISMOC, MSDS, IRIS-DMC and PDS for providing SEED SEIS data. F.K., N.S., Q.H. and D.K. were supported by NASA grant #80NSSC18K1628. V.L. acknowledges support by the Packard Foundation. C.N. was supported by strategic funds from the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We thank two anonymous reviewers for their constructive comments. This is InSight contribution 182.

The spectral envelopes of example events for each of four different event types defined according to their frequency content. The envelopes are obtained through visual selection on the event spectrograms. The event envelope is shown in black, whereas the S wave arrival and coda decay that is used in our analysis are shown in red. The gray thin line corresponds to a part of the data that is not used in our analysis, corresponding either to the noise level or glitches and features that are not part of the event signal.

Figure 2

The complete dataset of seismic events used in this study for vertical velocity seismograms (gray)

and envelopes (thick line). The seismograms and the respective envelopes are organized by the epicentral distance of the events. The black color represents the time window of the selected signal, whereas the part of the S-coda decay that was used in our analysis is shown in red.

Figure 3

Schematic showing the approximate average seismic raypath from the source to the station. The seismic waves that are produced at the source are scattered in the diffusive layer due to the presence of scatterers. These "real" ray paths are represented by the white rays, which are scattered in the diffusive layer. The approximation of the ensemble of these rays is represented by the blue curve, which is not scattered, but corresponds to a lower-than-true, apparent velocity in the diffusive layer.

Q = 400 are given for a range of layer thickness h = 1 -60 km with each color corresponding to a different V d /V e velocity ratio.

Figure 13

The envelopes obtained with the analytical method in this study (red, orange) are compared to envelopes computed through 2D numerical wave propagation modeling (black) [START_REF] Van Driel | High frequency seismic events on mars observed by insight[END_REF]. In the left panel, we use the same parameters and a range for the diffusivity, D = 0.5 -0.7 km 2 /s (i.e. the estimate of van Driel et al. ( 2021)) to compute the envelopes. The S-coda decay part of the envelopes show only a partial fit. In the right panel, we widen the range of diffusivities used to D = 0.3 -0.7 km 2 /s and find much-improved fits. The corresponding diffusivities are shown and we observe an increasing diffusivity that satisfies the data envelope as the ray path travels a longer distance in the diffusive layer. c) The results for Q = 400 are given for a range of layer thickness h = 1 -60 km with each color corresponding to a different V d /V e velocity ratio.